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CHAPTER  I 


INTRODUCTION 


This  report  describes  modeling  of  imploding  plasma  radi¬ 
ation  sources  (at  various  levels  of  description)  during  the 
contract  year  1982.  The  report  is  divided  into  two  parts. 
Part  I  describes  zero-D  models,  which  use  ordinary  differen¬ 
tial  equations  (and  their  numerical  or  approximate  analyti¬ 
cal  solution)  to  describe  the  dynamics  and  radiation  from 
imploding  wire  or  gas-puff  plasma,  with  the  approximation 
that  the  plasmas  are  uniform.  For  completeness,  Part  I  also 
includes  a  preliminary  description  of  some  results  actually 
obtained  at  the  beginning  of  the  1983  contract  year.1 
These  refer  to  first-order  nonadiabatic  corrections  to  for¬ 
mulas  for  radiation  from  quasi-adiabatic  pinch  plasmas. 

Part  I  concludes  with  a  report  on  analysis  of  a  ID  ide¬ 
alized  model  for  nonlinear  field  penetration  into  profile- 
invariant  plasma  where  the  conductivity  (i.e.,  field  diffu¬ 
sion  rate)  depends  on  field  strength  through  Joule  heating, 
etc  . 

Part  II  describes  ID  MHD  modeling  (or  more  exactly,  EMHD 
modeling,  since  the  electric  field  E,  rather  than  magnetic 
field,  is  explicitly  advanced  in  the  code),  with  related 
discoveries  about  the  nonlinear  physics  and  the  numerical 
techniques  for  advancing  the  variables.  Much  of  this  code 


1.  The  contract  reported  on  here  and  NRL  Contract  N00173- 
80-C-0202  have  been  merged  into  one  new  contract, 
DNA001 -83-C-0204  covering  02-03-83  to  03-31-84. 


development  has  been  carried  out  so  as  to  be  directly  usable 
in  2D  as  well,  although  the  2D  r-z  code  is  not  yet  near 
completion . 

The  purpose  of  the  formulations  described  in  Part  I  is 
to  provide  semi -quant i tat i ve  insight  and  understanding  into 
the  behavior  and  scaling  of  plasma  radiation  source  behav¬ 
ior,  using  simplified  model  assumptions.  The  purpose  of  the 
code  and  physics  described  in  Part  II  is  to  compute  quanti¬ 
tatively  the  behavior  of  plasma  radiation  sources  with  ade¬ 
quate  cylindrical  symmetry  and  to  provide  the  refinements 
and  modifications  to  the  simpler  less  quantitative  scaling 
laws  discussed  in  Part  I,  including  now  (e.g.,)  the  radial 
profiles  and  more  detailed  radiation  transport  which  are 
omitted  in  Part  I  for  simplicity. 

A  simple  "scoping"  code  is  discussed  in  Chapter  II  of 
Part  I.  In  this  code  called  RUNIN,  the  plasma  motion,  heat¬ 
ing  and  current  are  described  by  a  set  of  ordinary  differen¬ 
tial  equations.  The  code  follows  three  (or  two)  successive 
stages  in  the  evolution  of  a  wire  array  (or  annular  gas 
puff)  discharge:  the  motion  and  expansion  of  individual 
wire  plasmas  until  coalescence  to  an  annulus,  the  motion  of 
the  annular  plasma  (with  finite,  time-dependent  thickness, 
radius  and  temperature)  until  formation  of  a  non-hollow 
pinch,  and  the  compression  of  the  assembled  pinch  once  the 
annulus  closes.  In  each  stage  the  model  plasma  has  uniform 
density  and  temperature.  Radiation  from  the  wire  plasmas  is 
assumed  to  be  blackbody.  The  annular  and  pinch  phases  use 
blackbody-limited  otherwise-transparent  radiative  cooling 
rates  based  on  col 1 i s i ona 1-rad iat i ve  equilibrium  calcula¬ 
tions  without  doing  radiative  transport.  The  physical  basis 
of  the  code  was  discussed  in  JAYCOR  Report  J207-82-009  . 
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Simple  scaling  laws  are  derived  in  Chapter  III  of  Part  I 
for  the  radiation  from  a  quas i -ad iaba t ically  bouncing  pinch, 
in  which  the  radiative  loss  is  a  small  fraction  of  the  peak 
thermal  energy.  The  dynamics  of  the  wire-plasma  and  annular 
phases  are  also  discussed,  with  some  formulas  given  to  de¬ 
scribe  initial  conditions  for  the  pinch  phase.  The  scaling 
laws  are  used  on  uniform,  opt  ically-thin  pinch  plasma  mod¬ 
els,  and  corroborate  I"  dependence  of  the  energetic  radiated 
energy  (as  observed  at  Physics  International)  over  a  limited 
range  of  currents  and  plasma  masses,  (I«current  through  the 
plasma).  A  first-order  correction  to  the  minimum  radius  and 
maximum  temperature  of  the  implosion  is  also  given,  for 
small  but  non-negl ig i ble  radiative  energy  loss  dominated  by 
energetic,  transparent  radiation.  Caveats  to  the  model  are 
described . 

A  simple  model  for  the  nonlinear  penetration  of  electric 
field  into  a  finite  plasma  is  discussed  in  Chapter  IV  of 
Part  I,  in  which  the  conductivity  (and  thus  the  diffusion 
coefficient)  is  altered  locally  in  time  and  space  by  heating 
and  ionization  induced  by  the  diffusing  field.  A  ID  comput¬ 
er  simulation  is  described,  but  it  is  concluded  that  the 
restrictions  on  the  model  are  too  severe  to  accurately  por¬ 
tray  the  realities  of  field  penetration  in  imploding  annular 
plasmas,  and  the  reader  is  referred  to  the  more  advanced 
computational  results  of  Part  II. 

This  effort  is  intended  to  be  primarily  in  support  of 
the  DMA  objectives  of  understanding  and  projecting  an  ap¬ 
proximate  scaling  of  plasma  radiation  sources.  An  explicit 
tie-in  with  radiation  physics  work  funded  by  DNA  at  the 
Naval  Research  Laboratory  (NRL)  can  be  seen  from  the  central 


role  played  in  this  analysis  by  the  curve  of  energetic-pho¬ 
ton  radiative  power  as  a  function  of  plasma  electron  temper¬ 
ature  (it  also  depends  weakly  on  plasma  density  and  photon 
density)  -  such  curves  are  supplied  by  the  NRL  effort.  Les3 
visible  in  this  simplified  analysis,  but  also  important  in 
the  more  quantitative  calculations  enabled  by  Part  II  of 
this  report,  are  the  effects  of  radiation  transport  on  the 
energetics  and  the  radial  profiles  of  temperature,  etc.,  in 
imploding  plasmas.  Here  again  there  is  interaction  with  the 
plasma  radiation  physics  work  at  NRL,  both  in  the  use  of 
transport  algorithms  developed  at  NRL,  and  in  supplying  NRL 
with  reasonable  implosion  density  and  temperature  histories 
for  post -process ing  of  emission  spectral  histories,  etc. 
The  RUNIN  scoping  code  described  in  Part  I  of  this  report 
has  been  supplied  to  NRL,  to  use  in  conjunction  with  its 
more  extensive  radiation  physics  codes  as  an  Inexpensive, 
approximate  time-evolution  package  for  densities  and  temper¬ 
atures  during  implosions. 

Formulas  are  given  for  energetic  radiative  yield  (Y>)  as 

dependent  on  load  mass  (m)  and  length  (1),  experimentally 

measurable  current  in  the  load  (I),  kinetic  energy  at  the 

beginning  of  the  assembled  pinch  phase  (e,),  and  atomic 

physics  parameters,  such  as  the  temperature  (T  )  at  which 

P  K 

the  energetic  radiation  is  maximized. 

By  choosing  either  the  load  mass  per  unit  length  or  the 
current  so  as  to  optimize  Y>,  one  derives  a  formula  for 
scaling  of  optimal  energetic  yield  with  load  current  I,  in 
the  adiabatic  limit  where  the  radiative  energy  loss  has  only 
a  very  small  effect  on  the  dynamics.  Subject  to  the  caveats 
described  in  the  report,  we  verify  the  scaling 


found  by  earlier,  more  empirical  work  at  Physics  Interna- 
t ional . 

We  also  show,  to  first  order,  how  this  scaling  of  Y> 
were  I  falls  below  I"  at  higher  currents,  due  to  the  effects 
of  radiation  in  reducing  the  peak  attainable  temperature  for 
any  given  pinch  initial  conditions.  Here  ’first  order’ 
refers  to  the  quantity  Y ( rffl )/( thermal  energy),  Y(rm)  being 
the  total  radiative  energy  (energetic  +  less  energetic  pho¬ 
tons)  lost  between  onset  of  the  assembled  pinch  phase  and 
the  time  of  minimum  radius  r  (i.e.,  peak  compression). 


CHAPTER  II 


DEVELOPMENT  OF  THE  RUNIN  CODE 


For  the  assessment  of  Plasma  Radiation  Source  design 
feasibility,  it  is  important  to  have  an  approximate  scoping 
model  that  predicts  implosion  histories  and  scaling  of  radi¬ 
ative  output  with  the  appropriate  input  design  parameters. 
To  this  end  the  RUNIN  code  has  been  written,  incorporat  ing 
(with  modifications)  the  SAI  "WIRES"  Code  describing  the 
dynamics  of  an  imploding  array  of  wires,  the  JAYCOR 
"SQUEEZE"  Code  describing  the  dynamics  and  energetics  of  a 
pinching  cylindrical  plasma  (with  radiation  the  CRE  ioniza¬ 
tion),  and  a  "BRIDGE"  code  describing  the  intermediate  annu¬ 
lar  plasma  stage.  For  use  in  gas  puff  plasma  modeling,  the 
code  can  be  started  in  the  BRIDGE  phase.  This  RUNIN  code 
has  been  written,  documented,  run  and  debugged.  It  can  also 
serve  as  the  nucleus  of  a  variational  design  optimization 
code  1  . 

A.  WIRES:  EQUATIONS  FOR  THE  IMPLOSION  OF  THE  WIRE  ARRAY  - 

MODIFICATIONS  AND  UNDERLYING  ASSUMPTIONS 

1.  The  equations  governing  the  wire-plasma  radius  a 
(its  thickness,  not  its  position)  in  the  SAI  WIRES  code  was 
equivalent  to  the  assumption  that  the  radius  has  that  value 
which  gives  pressure  balance  between  the  thermal  pressure  of 
the  wire  and  the  magnetic  pressure  of  its  field.  The  tem¬ 
perature  is  determined  by  assuming  that  the  ZitaUc-T*  black- 

o 

body  radiative  cooling  exactly  balances  the  IZR  »  nJz£ira2 


8 


■US/. 


ohmic  heating  (baaed  on  full  field  penetration  of  the 
'wire').  Here  oB  is  the  Stefan-Boltzman  constant,  2ira2,  is 

D 

the  surface  area  of  a  wire  carrying  current  1^  over  cross- 
sectional  area  ira1  with  resistivity  n  (resistance  R).  We 
have  demonstrated  that  this  gives  instantaneously  quite  fat 
'wires'  at  early  times  if  I  is  large  and  T  small.  Since  the 
expansion  rate  of  the  wires  is  limited  by  one  or  a  few 
times  the  instantaneous  sound  speed,  we  have  limited  the 
radius  to  a„  ♦  / c  (t)dt  when  this  is  less  than  the  pressure- 
balance  radius.  Besides  being  more  physically  correct,  this 
prevents  the  code  from  merging  the  wires  into  an  annular 
plasma  at  very  early  times  before  they  could  have  merged 
physically . 

2.  When  the  N  wires  with  centers  at  radius  r  do  fi- 

c 

nally  merge,  i.e.,  when 

rcsin  ^  -  a  [II .  1 ] 

(the  tangency  condition  for  a  circular  array  of  N  circles  of 
radius  a),  they  are  replaced  by  an  annular  layer  (with  outer 
radius  rt  and  inner  radius  r,)  which  has  the  same  area  as 
the  collection  of  wires 

ir(r?  -  r?)  -  Nira*  [II  .2] 

and  has  a  geometric  mean  radius  equal  to  the  radius  of  cen¬ 
ters  of  the  wires:  r.r,  -  r  2.  The  new  velocity  of  the 

'  c 

outer  surface  rlf  is  given  by  conservation  of  kinetic  energy 
in  the  transition. 

The  use  of  the  notation  r,  for  the  outer  radius  and  r5 
for  the  inner  radius  of  the  resulting  annulus  is  intended 
to  conform  (at  least  in  subscripts)  to  the  notation  of  the 
RUNIN  code,  in  which  y,  and  y.  are  the  outer  and  inner  radii 


respectively.  The  code  announces  the  merge  (both  on  the 
terminal  and  in  the  output  data  file)  and  proceeds  to  the 
integration  of  equations  of  motion  for  a  uniform  annular 
plasma,  using  a  high-precision  package  integration  routine. 

The  magnetic  field  energy  in  the  WIRES  code  includes  a 
contribution  from  positions  inside  the  wire  array.  This  is 
physical,  but  should  go  to  zero  as  the  wires  merge  to  an 
annulus,  even  if  the  approximate  equations  of  the  WIRES  code 
break  down  when  the  wires  are  close  together.  In  the  RUNIN 
code  we  have  taken  this  contribution  zero  as  soon  as  the 
wires  merge. 

Consistent  with  the  assumption  of  complete  field  pene¬ 
tration  of  the  individual  wires,  we  have  assumed  uniform 
current  density  in  the  annular  layer  when  it  comes  to  calcu¬ 
lating  energy  density  integrals,  but  we  have  assumed  in  the 
force  equations  that  all  the  J  x  B  force  is  applied  to  the 
outer  surface.  In  fact,  a  real  physical  difference  between 
wire  array  implosions  and  foil  or  annular  gas  puff  implo¬ 
sions  may  be  that  for  wire  array  implosions  the  field  pene¬ 
tration  is  fairly  complete  early  in  the  wires  phase  and 
continues  to  be  fairly  complete  after  merging  to  an  annulus, 
while  in  foil  or  gas  puff  implosions  the  current  may  be 
restricted  to  a  relatively  thin  layer  near  the  outer  surface 
of  the  annulus.  The  full  ID  MHD  code  described  [elsewhere 
in  this  report]  is  capable  of  estimating  the  degree  of  field 
penetration  into  foils  or  gas  puff  annuli,  and  of  showing 
the  sensitivity  of  the  implosion  dynamics  to  the  degree  of 
penetration,  whereas  the  simpler  model  presented  in  this 
section  cannot  investigate  that  question  consistently.  The 
actual  field  penetration  is  discussed  in  more  detail  in 
Chapter  II  of  Part  II  of  this  report. 


B.  BRIDGE:  EQUATIONS  FOR  THE  UNIFORM  ANNULAR  PLASMA  IMPLO¬ 
SION 

The  equations  of  motion  of  the  annular  plasma  result 
from  performing  /dr  on  the  fluid  equation  and  assuming  that 
the  density  gradient  3n/3r  remains  zero  except  at  the  plasma 
edges.  This  latter  assumption,  equivalent  to  3/3r(V«v)  -  0, 
requires  v(r)  of  the  form 


v(r )  -  k ,r  +  kar 


[II. 31 


and  the  conditions  v(r,)  =  rlt  v(r,)=rs  give  equations  for 

kj(rlt  r,,  rs,  rs)  and  k,(r, ,  rt,  rs,  r5): 


1  r  t 

klp1  *  W*  r1  -  W7  r5 


k2/r5  "  “  W7  r1  +  TrrT  r5* 


[II. 4] 


[II. 5] 


Here  r=r5/rlf  the  ratio  of  inner  to  outer  radius,  which  is 

of  course  time-dependent  (as  are  kv  and  ka).  In  the  inte- 

2 

grated  force  equation  (Tidman  and  Colombant,  1979)  one  has 


then 


/r'  7T  dr  ♦  /r‘  vVvdr  -  \  kt(rf  -  r  , ) +k ’ In ( r , /r s ) 
J  r  s  9 1  J  r ,  2  1  1  *  2  11 


♦  i  (r?  -  rf) 


[II. 6] 


The  internal  viscous  stress  tensor  the  the  velocity 
field  of  Equation  [II. 3]  is  divergenceless,  so  that  all  the 
force  comes  from  the  thermal  and  magnetic  pressure  gradi¬ 


ents.  For  the  net  force  we  have 


/'  1  VPdr  -  (n  ♦  n.  )  T 
r  j  e  l 


B 2  ( r  ,  ) 


[II. 7] 


with  B(r, )  -  2I/cr ,  (in  cgs  units). 


r  j  f  V  .*  «*  /.V.  .•  ■  . 


It  will  be  seen  that  if  the  velocity  of  the  inner  sur¬ 


face  remains  finite  as  it  collapses  to  the  origin,  then  the 
in  term  in  Equation  [II. 6]  is  singular  as  rs  -*•  0.  It  is 
thus  a  mathematical  feature  of  these  approximate  equations 
that  a  finite  excess  pressure  from  the  outside  causes  rs  to 
take  on  a  momentarily  infinite  negative  value  as  r,  •+•  0. 
This  lasts  for  zero  duration  and  causes  no  physiol  discrep¬ 
ancy  because  the  singularity  is  integrable  and  rt  remains 
finite.  It  does,  however,  mean  that  the  numerical  integra¬ 
tion  of  the  equations  must  be  done  cleverly  and  approxi¬ 
mately  near  the  time  of  transition  from  annulus  to  plasma 
cylinder.  Since  the  integrating  routine  used  (DGEAR) 
shrinks  its  timestep  when  it  finds  larger  accelerations, 
|r,|  is  arbitrarily  limited  to  a  reasonable  value  (10,7cm/ 
sec2)  and  the  equations  are  carried  forward  approximately 
for  the  remaining  short  (sub-nanosecond)  time  until  rs  •  0. 

C.  NUMERICAL  INTEGRATOR  PACKAGE 

While  the  WIRES  subroutine  calls  a  simple  fixed-step 
integrator  subroutine  (the  3ame  as  in  the  original  SAI  WIRES 
code),  the  BRIDGE  and  SQUEEZE  subroutines  use  the  IMSL  Li¬ 
brary  DGEAR  package,  a  var iable-timestep  pred ictor-corrector 
integration  routine  designed  for  use  with  sets  of  stiff 

3 

differential  equations  .  An  override  and  exit  is  provided 
when  doinjg  the  integrations  in  the  BRIDGE  phase,  so  that  the 
inner  radius  toes  to  zero  within  some  tolerance  (10~*cm) 
without  going  negative  or  shrinking  the  timestep 
indefinitely . 


a  /•;%  .**  -  .*•  •  -  .*«.  ■  .  -  •  .  • 
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iid 


i/2irrdr  nf(r)  Z{T(p)}  F_  „  „  { T  ( r  ) } 


in  the  optically  thin  limit, 


pad  '  2*^Veff 
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in  the  blackbody  limit,  and  not  exppessible  in  closed  fopm 
except  in  these  limits. 

Yet  a  simplified  implosion  model,  to  be  useful  and  show 
coppect  trends,  must  include  in  some  approximate  way  an 
average  radiation  function  with  about  the  same  magnitude  and 
dependence  as  the  more  exact  problems,  i.e.,  a  suitable  set 
of  parameters  n^O),  rR ,  T(0),  r^,  T(«)  must  be  chosen, 
representing  the  heights  and  widths  for  the  n^r)  and  T(r) 
distributions,  and  a  function  Rad(r,,  rt,  T)  must  be  con¬ 
structed,  so  that 


irUCrJ  -  r\)  n*Z(T)Rad(r , ,  rs,  T)  -  P( 


[11.11] 


for  suitable  n^n^O),  r^l  ,  T{T(0),  T(«),  i"T)  •  The  proper 
choice  of  such  parameters  and  suitable  distributions  for 
density  and  temperature  are  being  investigated  using  the 
BLANKET  code. 

For  the  radiative  output  above  1  keV  the  effect  of  non- 
isothermal  temperature  profiles  in  the  assembling  plasma  is 
indicated  semi  quant i tat i vely  by  using  the  transparent  CRE 
model  with  Gauss  ian-plus-constant  temperature  profiles  in  a 
Gaussian  density  profile.  Figure  II. 1  shows  the  ratio  of 
radiated  power  above  1  keV  for  this  distribution,  in  alumi¬ 
num,  to  that  for  an  isothermal  aluminum  plasma  with  the  same 
Gaussian  density  profile,  as  a  function  of  the  relative 
width  (r^)  of  the  temperature  profile.  The  temperature 
distributions  are  chosen  in  this  case  to  have  the  same  cen¬ 
tral  temperature,  and  the  temperature  of  the  cooler  outer 
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Ratio  of  radiated  power  above  1  keV  for  Gaussian-plus-constant 
temperature  distribution  to  that  from  a  uniform  temperature 
distribution.  Both  distributions  are  assumed  for  plasma  cyl¬ 
inder  with  Gaussian  density  profiles  of  unit  width,  and  the 
parameter  r measures  the  relative  width  of  the  temperature 
distribution  =  «  for  isothermal  limit) . 


T  1  s  •  l  ■  r  •  i  ■  ■  1  i  i  1  ■  .■  i'  .  r  j  — 'j-i-i 


region  is  varied  as  a  parameter.  For  use  with  the  RUNIN 
code,  where  the  energetics  are  described  by  a  single  temper¬ 
ature,  comparisons  are  made  with  the  same  average  tempera¬ 
ture  rather  than  with  the  same  central  temperature. 

E.  CAVEATS  TO  THE  RUNIN  CODE 

The  code  described  above  models  the  dynamics  of  uniform 
plasma,  either  as  (1)  identical  symmetrically-placed  wire 
plasmas  of  c lrcular  cross  section  with  forces  acting  at  ! 

their  centers,  or  (2)  as  an  annulus  of  uniform  density  with  j 

force  acting  on  the  outer  surface  and  model  equations  for  * 

the  motion  of  the  inner  radius,  or  (3)  a3  a  uniform  z-pinch 
driven  at  its  outer  surface  by  the  net  pressure  (thermal  t 

minus  magnetic)  but  ohmically  heated  as  though  carrying  ! 

uniform  current  (a  feature  easily  generalized  if  the  cur-  * 

rent-carry ing  area  is  known).  The  circuit  equation  coupling  ! 

the  plasma  load  current  to  the  generator  voltage  includes 
diode  inductance  and  resistive  impedance  but  does  not  pre¬ 
sently  include  any  currents  not  drawn  in  the  load.  The 
radiative  loss  is  assumed  either  transparent  or  (when  this 
exceeds  blackbody)  blackbody-limited .  In  the  pinch  phase 
this  is  probably  a  fair  description,  although  the  assump¬ 
tion  of  uniform  temperature  and  density  gives  an  overesti¬ 
mate  of  the  energetic  radiation  (as  just  discussed  in  Sec¬ 
tion  D).  In  the  case  of  annular  plasmas  with  moderately 
high  mass  and  large  radius  the  radiation  formula  is  less 
believable  and  tends  to  over-cool  the  plasma. 


Time-integrated  pinhole  and  streak  photos  of  wire  array 
implosions  appear  to  indicate  that  for  massive  enough  arrays 
there  is  considerable  plasma  blowoff  inward  from  the  denser 


wire  plasmas,  forming  either  low-density  plasma  implosions 
or  plumes^.  This  phenomenon,  though  interesting,  is  not 
within  the  scope  of  the  WIRES  model,  which  assumes  uniform 
wire  plasmas.  Because  of  the  observed  advantages  of  gas- 
puff  implosions,  this  wire  blowoff  phenomenon  has  not  been 
further  investigated  here,  even  though  it  represents  a  sig¬ 
nificant  deviation  from  the  WIRES  model. 

One-D  MHD  codes  allowing  radial  resolution  sometimes 

show  the  formation  of  a  density  and  temperature  spike  on  the 
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axis  when  the  plasma  annulus  closes  .  This  is  entirely 
reasonable  and  physical,  but  represents  another  feature  not 
describable  by  a  uniform-plasma  model.  In  this  case,  how¬ 
ever,  we  hope  to  generalize  the  simplified  model  to  allow 
such  a  central  compressed  zone  because  the  dynamics  of  gas 
puff  runs  may  be  seriously  affected  a3  the  dense  central 
spike  is  radiatively  cooled  to  allow  further  compression. 


CHAPTER  III 


APPROXIMATE  SCALING  LAWS 


A.  WIRE  PLASMA  PHASE  OF  ARRAY  IMPLOSIONS 


The  very-early-time  behavior  of  wire  array  implosions, 
involving  the  phase  transitions  of  the  material  and  the 
associated  "pause"  (dip)  in  the  conductivity,  is  discussed 
in  the  literature  on  exploding  wires  and  is  not  treated 
here.  For  practical  modeling  of  the  implosion  dynamics,  we 
start  with  an  initial  temperature  at  which  the  plasma  has  no 
more  neutrals  and  is  not  dominated  by  strong  coupling  but 
has  classical  conductivity  determined  by  Coulomb  collisions. 

After  an  initial  heating  expansion,  the  individual  wire 
plasmas  reach  a  local  dynamic  thermal  equilibrium  in  their 
own  reference  frames.  For  the  case  of  N  imploding  wires, 
local  Bennett  equilibrium,  (i.e.,  pressure  balance)  gives  an 
equilibrium  temperature  T: 

m.  I  * ( tSt  ,  )  m 

I1*Z(T)]I  .  •  ji  6(t)  CIII.U 
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where  m^  is  the  mass  of  an  ion,  v  is  the  total  load  mass 
(all  N  wires)  per  unit  length,  Affl  is  the  atomic  number  of 
the  load  and  IMA  is  the  total  current  (in  MA)  in  all  the 
wires.  The  degree  of  ionization  Z(T)  for  aluminum  at  densi¬ 
ties  of  order  1021  ions/cm*  is  shown  in  Figure  III.l.  The 


18 


individual  wire-plasma  radius  p  expands  until  blackbody 
radiative  loss  equals  the  ohmic  heating  power  and  prevents 
further  heating  except  by  increasing  current.  This  occurs 
at  a  wire-plasma  radius  given  by 


C2.pl)(oBT*>  -  l^)  (i)‘ 


[III. 3] 


with  o  the  conductivity.  Once  the  temperature  is  high 


enough  that  o  is  determined  by  Coulomb  collisions. 
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(so  defining  the  constant  K  ) .  When  these  wire  plasmas  just 
merge,  their  centers  are  at  radius  pcsc(ir/N)  and  the  result¬ 
ing  annular  layer  has  outer  radius  given  by 


r\  -  [(M/2)  +  /(N/2)1  +  CSC w  ( n/N )  ]  p? 
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with 
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in  terms  of  the  current  I,  and  temperature  T,  at  time  tt 
(i.e.,  merge).  (The  geometrical  factor  in  brackets  is  8  for 
6  wires  and  22  for  12  wires.)  Numerically, 


,__x  ,,o  (ll  (ma)  ^2/3,,„  si/3„-n/6 

i  (  cm )  -  172  l — jj - J  Z(T,)  Ti(ev)> 
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and  using  the  preceding  information  for  T ( I )  in  this  expres¬ 
sion  one  ha3 

{  1  >Z  (  T  ,  )  }  l11/67/6 

P  »  (  c ® )  *  .123  I » (  M A  )  A  ^  {  Z  (  T  ,  )  ) 
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But  this  scaling  of  p4  and  r,  with  I-1  [or  about  I-1’6  when 
the  Z ( T ( I ) )  dependence  is  included]  holds  only  over  a  limit¬ 
ed  range  of  load  mass  and  current. 

As  long  as  the  wire-plasma  radius  p  is  not  clamped  by 
blackbody  loss,  it  expands  at  a  rate  proportional  to  the 
sound  speed,  i.e.,  p  =  {(1+Z(T)}1/2  I  «  In1'2,  until  the 
plasmas  merge.  With  p  <*Ip±l/2  and  the  radius  of  centers  rQ 
moving  as  rcr'c  I  2  /  y ,  the  merge  condition  p»r  c  sin(ir/N) 
determines  the  annulus  formation  in  a  dynamic  way,  if  this 
merge  occurs  before  radiation  can  limit  p.  For  constant 
current  the  merge  condition  for  non-equ i 1 i brated  p  (i.e., 
for  high  current)  is 


.  .  P  -  /771  r0/2N 
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where  r0  is  the  initial  array  radius.  For  N«6  wires  this 
gives  p-0.1  r0  and  rc  *  0.3  r  „  at  merge.  Thus,  if  radia¬ 
tion  does  not  limit  the  wire-plasma  radii,  the  merge  occurs 
before  very  significant  implosion  has  taken  place,  because 
the  superheated  wires  expand  rapidly.  One  then  has  an  annu¬ 
lar  plasma  implosion  (starting  with  inner  radius/outer  radi¬ 
us  =  1/2  for  6  wires)  for  more  of  the  implosion  than  if  the 
wire  plasma  size  had  been  radiation  limited. 

From  the  code  and  the  scaling  just  discussed,  we  can 
identify  the  following  partially  overlapping  stages  in  the 
WIRES  phase  of  a  wire-array  implosion: 


(1)  Ohmic  heating  dominates;  wires  expand  and  heat; 


(2)  Ohmic  heating  ■  PdV  cooling  (a  p);  wires  expand; 


(3)  Bennett  equilibrium;  P„  -P  >0;  (1+Z(T)}T  *  I2/yN.  Radi- 

U  U  D 

ative  cooling  adjusts  wire  radius  if  T  is  low  enough  for 
the  radiation  to  be  blackbody; 

(4)  Blackbody  /  Bennett  equilibrium:  {1*Z(T)}T  31  I2/yN  and 

„  „  7,t.1/3t2/3t-11/6m7/6  1  1  /  6  _ -  1  .6  7/6 

p  11  Z  (  T  )  I  T  N  ,  u  I  N 

(Not  applicable  if  T  becomes  high  enough  for  quasi¬ 
transparent  radiative  loss.) 

Merging  of  the  wire  plasmas  may  interrupt  the  process  at 
stage  (2),  (3)  or  (4),  or  there  may  be  a  further  stage  if 
the  wires  overheat. 

In  the  simple  physics  of  the  WIRES  phase  of  the  RUNIN 
code,  only  stages  (1)  and  (4)  are  explicitly  included. 
Figure  III. 2  shows  these  stages,  occupying  roughly  equal 
fractions  of  the  WIRES  phase  at  constant  driving  voltage. 

The  modeling  and  scaling  of  this  wires  phase  is  impor¬ 
tant  to  the  radiation  and  dynamics  in  the  later  pinch  phase 
only  through  the  way  it  scales  the  initial  conditions  of  the 
pinch  phase  (and  in  general  through  the  duration  of  the 
wires  phase,  if  the  generator  pulse  time  is  limited). 

The  radiative  yield  of  the  pinch  depends  on  the  initial 
radius  and  temperature  (at  the  beginning  of  the  pinch  phase) 
as  r_l  T  ~ 1  /  2  .  While  I2/y  is  low  enough,  this  goes  as  y_,,/3 
I2. 

B.  ANNULAR  PLASMA  STAGE 

The  merging  of  wires  produces  a  somewhat  irregular  annu¬ 
lar  plasma;  ignoring  the  i rregular  it  ies ,  one  can  derive 
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Figure  III. 2.  Radial  positions ,  current  and  energetic  line  radiation  versus 
time  (in  ns)  for  a  typical  RUNIN  code  run.  The  wires  expand 
until  radiative  cooling  balances  ohmic  heating  (dashed  verti¬ 
cal  line  marks  the  onset  of  such  Blackbody-Bennett  equilib¬ 
rium)  .  Inner  and  outer  extent  of  the  wires  are  shown  in  the 
WIRES  stage  (radius  of  centers  shown  dashed);  inner  and  outer 
radii  of  the  annulus  are  shown  in  the  BRIDGE  stage.  Pinch 
stage  (SQUEEZE)  onset  occurs  when  inner  radius  goes  to  zero. 
The  run  shown  was  for  constant  3  MV  voltage . 
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differential  equations  for  the  inner  and  outer  radii  and  the 
temperature  in  various  simplified  models.  The  same  models 
and  equations  can  be  applied  to  the  annular  plasma  formed  by 
the  pinching  of  an  annular  gas  puff  from  a  supersonic  ring 
nozzle.  Even  with  the  simplest  physical  models,  however,  it 
has  not  been  possible  to  do  the  sort  of  analytic  modeling 
here  that  has  been  possible  for  the  later  pinch  phase  of  the 
implosion.  Here  then  the  RUNIN  code,  or  specifically  its 
BRIDGE  subroutine,  provides  the  only  simple  guide.  The  code 
has  been  used  to  show,  for  example,  that  the  ratio  of  final 
to  initial  radii  and  temperatures  for  the  annular  phase  are 
very  insensitive  to  the  initial  radius. 

To  develop  simple  rules  for  expressing  variables  at  the 
end  of  the  annular  phase  in  terms  of  those  at  the  beginning 
of  the  annular  phase,  the  BRIDGE  subroutine  was  run  for  a 
range  of  input  variables.  The  "output  variables"  chosen  for 
examination  were  ratios  of  final  to  initial  radius,  tempera¬ 
ture,  kinetic  energy,  and  r/T,  a  variable  which  occurs  in 
the  formula  for  radiative  energy  loss  in  the  case  of  trans¬ 
parent  radiation. 

Table  III.1  shows  the  percentage  variation  in  the  output 
variables  divided  by  the  dynamic  range  (maximum  value/ 
minimum  value)  of  the  input  variables,  in  the  absence  of 
radiation.  None  of  the  output  variables  depend  sensitively 
on  initial  radius  r,.  The  kinetic  energy  ratio  is  the 
only  output  variable  to  have  an  interesting  dependence  on 
the  initial  velocity  v  of  the  annulus  outer  radius  (i.e., 
on  the  initial  kinetic  energy  of  the  annulus);  the  ratio 
of  initial  to  final  values  of  r,  T,  and  especially  r/T  are 
insensitive  to  vi  .  The  r/T  ratio  (final  to  initial,  for 


Dependent  Variables 


r .  Jt~. 
i  1 


0 

0 

0.6/2 

0 

-.16/30 

.13/30 

-318/30 

.09/30 

.18/4 

-.34/4 

-.67/4 

-.09/4 

Table  III.l. 


Variation  in  ratios  of  final  to  initial  dynamical  quantities 
(shown  by  numbers  in  numerators)  over  dynamic  range  of  ini¬ 
tial  conditions  for  the  plasma  annulus  (shown  by  factors  in 
denominators;  a  factor  of  30  in  v ^ ,  for  example,  indicates  a 
factor  of  30  change,  about  some  nominal  value  typical  of  the 
inward  speed  at  the  time  the  annulus  is  formed).  Subscript  f 
refers  to  the  time  at  which  the  center  of  the  annulus  closes. 
A  zero  indicates  very  weak  dependence. 


the  annular  phase)  is  also  insensitive  to  the  initial  tem¬ 
perature  or  the  annular  plasma,  ,  although  the  r  and  T 
ratios  depend  on  T^. 

C.  RADIATION  SCALING  AND  DYNAMICS  IN  THE  PINCH  PHASE 

1 .  The  Nearly-Adiabatic  Uniform  Radiator 

The  RUNIN  code  indicates  that  the  temperature  versus 
radius  behavior  of  the  pinch  phase  lies  fairly  close  to  the 
adiabatic  law,  for  low-ma33  load3  and  low  currents.  The 
rate  of  doing  PdV  work  greatly  exceeds  both  ohmic  heating 
and  radiative  loss  except  for  a  very  brief  interval  near 
minimum  radius.  In  addition  the  code  indicates,  in  agree¬ 
ment  with  experiments,  that  the  current  does  not  change  by 
a  large  factor  during  the  compression.  These  facts  suggest 
an  approximate  set  of  differential  equations  which  can  be 
treated  analytically,  with  some  approximations,  to  yield 
formulas  for  the  output  of  energetic  radiation  and  for  its 
optimization.  This  section  of  the  report  treats  first  the 
zero-order  formulas  that  result  for  a  pinch  in  the  adiaba¬ 
tic  limit,  and  then  the  first-order  correction  when  the 
radiative  energy  loss  deepens  the  implosion  perceptibly  but 
not  so  much  as  to  make  the  r-T  trajectory  deviate  far  from 
the  adiabatic  one  for  very  long  during  the  implosion.  The 
utility  of  these  formulas  lies  in  the  understanding  they 
provide  of  the  central  qualitative  concept  of  how  to  design 
and  optimize  an  imploding  radiator,  and  in  their  scaling 
with  several  variables  at  once  -  load  mass,  atomic  mass 
number,  initial  conditions  at  the  assembly  of  the  plasma 
pinch,  total  current  delivered  to  the  load,  etc.  They 


reflect  the  facts  that  overdriving  or  underdriving  the  load 


can  result  in  a  radiation  pulse  that  is  too  short  or  too 
cool,  or  occurs  at  a  less-than-optimal  density;  they  indi¬ 
cate  what  the  optimal  load  current  is  for  given  conditions 
at  the  formation  of  the  pinch,  and  how  it  depends  on  the 
load  variables,  how  sensitive  the  radiation  loss  is  to  devi¬ 
ation  from  optimum,  and  what  the  peak  attainable  yield  is  as 
a  function  of  the  variables,  subject  of  course  to  the  limi¬ 
tations  and  approx imat ions  of  the  model.  Finally,  a  discus¬ 
sion  is  provided  of  the  limitations  of  the  model,  and  of 
what  happens  to  the  predicted  versus  actually  expected  radi¬ 
ative  behavior. 

We  arrive  at  simple  formulas  in  two  ways:  first,  by 
simple  dimensional  evaluation  of  the  t ime- i ntegral  of  the 
radiated  power,  and  second,  by  a  more  careful  evaluation 
based  on  expansion  of  the  integrand.  Both  techniques  pre¬ 
suppose  that  the  plasma  is  nearly  transparent  to  the  ener¬ 
getic  radiation,  and  that  is  consistent  with  the  load  masses 
of  present  and  near-term  experiments.  We  note  that  similar 
techniques  apply  to  neutron  production  from  fusion  targets 
inasfar  as  the  targets  remain  uniform  cylindrical  or  spher¬ 
ical  plasmas,  but  to  radiatively  driven  implosions  only 
inasfar  as  the  radiation  pressure  varies  nearly  inversely 
with  load  radius. 

A  uniform  pinch  with  mass  m,  atomic  mass  m^ ,  radius  r, 
length  l,  and  temperature  T\-Te-T  is  governed  approximately 
by  the  dynamical  equation 


mr  -  2irr  £,(  P.  .  .  -  P  )  , 

thermal  mag 


with 


thermal 


n  £  ( 1 *Z)T, 


(m/m^  ) 
irr  *  4 


[III .  1  0] 
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[Ill . 1 2] 


Pmag  "  T:  [§F>'  in  083  units 

Here  Z(T)  is  the  degree  of  ionization  and  I  is  the  pinch 
current.  This  equation  has  the  form 

rr  -  aT  -  b  [II 1. 13] 


when  Z  and  I  are  nearly  unchanged  during  the  course  of  the 
pinch.  Here 


a  =  2(1 +Z)/a 


[III  - 1  4] 


and 


[III. 15] 


When  ohmic  heating  is  negligible  (as  with  short-duration 
pinches  with  high  conductivity)  the  temperature  is  governed 
approximately  by 

T  -  -{(Cr/r)  +  f } T  [III. 16] 


where  C  is  an  approximately  constant  adiabatic  index,  and 
the  radiative  loss  function  f(r,T)  will  be  neglected  in  the 
adiabatic  limit  and  used  as  a  perturbation  in  deriving  the 
first-order  corrections.  The  temperature  model  equation 
arises  from  equating  PthV  power,  less  radiative  loss,  with 
the  time  derivative  of  3 / 2 N ^ ( 1 *Z )T  +  Niei’  where  is  the 
mean  ionization  energy  (as  a  function  of  T)  and  N.-m/m..  is 
the  total  number  of  ions  being  compressed.  The  thermal 


while  V=2*irr.  When  f  i3  negligible  this  T  equation  can  be 
integrated  to  give  the  adiabatic  law, 


T/T,  -  (rt/r)  , 


[III .  1 8] 


with  r,  and  T,  the  initial  values  at  time  t,  when  the  pinch 
phase  begins.  This  form  for  T  is  substituted  into  the  force 
equation,  which  can  then  be  integrated  once  to  give 


Cr2  -  Cr 2  -  2a(T-T,)  +  2biln(T/Tl) 

or  its  equivalent  form  involving  r  instead  of  T. 
The  energetic  radiative  yield  is  the  integral 


[III .  1 9] 


(m/m.  )  2 


-7 


;{T(t)  }dt 
r  2  ( t ) 


[III  .20] 


in  the  transparent  limit.  The  function  g(T),  i.e.,  the 
energetic  radiation  per  pair  of  ions,  is  known  for  colli- 
sional  radiative  equilibrium,  and  it  is  assumed  that  the 
dynamical  changes  are  slow  enough  that  changes  in  radiative 
power  track  them  without  appreciable  delay.  The  irlr2  denom¬ 
inator  comes  from  the  assumption  of  uniform  density  and 
temperature,  and  gives  a  resonance  shape  to  the  integrand, 
which  peaks  near  the  time  of  peak  compression.  If  the  maxi¬ 
mum  temperature  T  is  near  the  peak,  T  .  ,  of  g(T),  then  the 

m  P* 


integrand  peaks  when 
minimum  radius. 


reaches 


and  r  reaches  r 


In  the  simplest  approximation  one  evaluates  the  integral 
using  the  resonance  approximation. 


f  SLilUt 


g(T  )x 
m 


[III  .21  ] 


where  the  pulse  width  t  is  found  approximately  by  treating 
r(t)  as  quadratic  about  peak  compression; 


1 


r  +  —  r  t 2 
m  2  m 


(r  /r  ) 
m  m 


1  /  2 


[III  .22] 


so  that 


T_ 

r 


r  / ( rr  ) 
m  m 


r  (aT  -b ) *  7  2 
m  m 


[III  .23] 


Generally  aT  >>b  during  most  of  the  pinch  phase.  To  the  ex- 

81  ■*  1  /  C 

tent  that  r(T)  follows  the  adiabatic  law,  one  has  r  «T  , 

m  m 


leaving  only  a  weak  dependence  of  t/r2  on  T  .  This  gives 


y>- 


(m/m^  ) 1 


irta1  /2r  ,T, 


-(1/C-1/2),  . 

T7C  Tm  8  Tm 


[III.2H] 


in  the  adiabatic  resonance-approximation  limit.  Recalling 


that  a-2(1+Z)/m.  and  using  Z>>1  and  Z/m,  »  1/2m  (with  m 
i  i  p  p 


the  proton  mass),  one  can  write 


.  u/n2  l».-2T,/C-,/2g(T  ) 
1  m  m 


[III  .25] 


To  the  extent  that  1/C-1/2  is  small,  this  is  clearly  maxi¬ 
mized  (for  given  rt,  Tt)  if  the  driving  current  and  initial 
conditions  are  such  as  to  make  the  maximum  temperature  Tm 
coincide  with  the  peak  of  the  energet ic-rad iat ion  curve 
g(T)  . 

Since  in  the  adiabatic  limit  T  is  reached  together  with 

m 

r^  when  r-0,  we  have  from  the  equation  for  Cr 2  an  equation 

for  T  : 
m 

a(Tm-T,)  -  b  tn ( T^/T ,  )  -  1  Crf,  [III. 26] 


A  useful  form  of  this  for  iteration  is 

Tm  "  Tl  +  n  +  a  ln(Ta/Tt>»  [III. 27] 

since  the  in  term  varies  slowly.  The  rf  term  represents 
kinetic  energy  provided  prior  to  t,  (which  may  also  depend 
on  a  and  b),  and  the  in  term  represents  conversion  of  mag¬ 
netic  field  energy  to  thermal  energy  since  time  t,. 

Prior  to  time  t,  the  motion  is  not  given  by  the  same 
differential  equations,  is  generally  not  adiabatic,  and 
generally  does  not  involve  near ly-constant  current.  The 
kinetic  energy  acquired  by  time  tt  does,  however,  scale  with 
b  (i.e.,  with  I*)  for  any  given  pulse-shape  of  I(t)  deliv¬ 
ered  to  the  load.  Thus  if  T  >>T, ,  then  T  is  approximately 
proportional  to  b,  i.e.,  to  t,I2/m.  So,  while  refinements 
can  and  will  be  made  to  this,  one  ha3  a  rough  rule  of  thumb 
that  for  given  r,  and  T,,  the  current  I  required  to  drive 
T  up  to  Tpk  scales  nearly  as  (yTpk)1'2,  where  u=m/l,  in  the 
adiabatic  limit;  and  Y>  scales  roughly  as  uzim~*g(  Tp|< )  . 

Both  scalings  are  equally  valid  if  I  .  is  replaced  by  I 

ope  m 

and  T  ,  by  T  ,  so  that  one  has 
pk  J  m 

Y>  «  I  ■*  im  “  2  g  (  T  )  T  ~2  [III. 28] 

i  mm 

(aside  from  In  terms),  in  the  adiabatic  limit.  With  less 

than  optimal  current,  i.e.,  where  T  is  on  the  shoulder  of 

m 

the  g(T)  curve  so  that  g(T)  is  increasing  roughly  as  a  power 

of  T,  one  may  even  find  a  region  where  g(T  )T-2  a  T°  and  so 

m  m  ^ 

experiments  in  this  domain  could  fall  along  a  curve  Y  * 
1*1,  even  at  constant  m. 

We  will  see  later  that  the  current  required  to  reach  a 

given  T  is  larger  in  the  nonadiabatic  case.  This  will  be 
m 

only  slightly  offset  by  the  fact  that  r  will  be  smaller, 


giving  larger  values  of  fg/r2  once  the  desired  T  is  reach- 

111 

ed.  In  the  totally  lossy  case,  the  plasma  cools  as  fast  as 
it  can  be  compressi vely  heated,  and  temperatures  on  the 
order  of  Tpk  cannot  be  reached.  In  that  case  Y>  is  general¬ 
ly  small.  This  can  happen  if  the  load  is  sufficiently 
massive . 

The  evaluation  and  optimization  of  Y >  can  be  refined 

somewhat  within  the  adiabatic  limit.  One  can  make  a  better 

approximation  to  /g(T)dt/r2  and  one  need  not  argue  that  T°pt 

-  T  .  To  evaluate  the  integral,  one  first  constructs  T/T 
P  K 

as  a  function  of  $-in(T/Tt),  using  the  expressions  for  r(T) 

and  r(T),  and  then  uses  this  to  change  the  integration 

variable  from  t  to  $.  One  then  expands  the  singular 

denominator  of  the  new  integrand  about  its  branch  point  at 

T-Tm,  using  the  method  of  fractional  derivatives  (rather 

than  simply  expanding  ra  as  r2  +  1/2  r  t2).  Then,  because 

m  m 

the  g ( T )  curve  is  nearly  parabolic  in  a  log ( g ) - log ( T )  plot, 
one  can  use  the  model 


ing  -  in  g  -  B[ in ( T/T  ,  )  ] 2 
p  k  pk 


[III. 29] 


(with  fit  parameter  B)  to  do  the  resulting  integral.  Car¬ 
rying  out  this  procedure,  we  have 


A* 


+  bi|>-aT , 


[III. 30] 


with  $  -  in(T/T,),  A  »  C~l  *  2B$p(<,  e,  -  ^  Crf,  and 

e,  ♦  b$  -  aT.Ce1*  -  1)  -  0  defining  $  .  [III. 31] 

in 


Letting  w-/$  -$  and  expanding  the  radicand  about  w-0, 


c  j  +  b$  -  aT  t  ( e 


b )  wz  , 


[III. 32] 


-  1 )  «  (aT  - 
m 

the  integral  on  the  right  becomes 


p<b,*a) 


/  (b ) 
m 

o 


Q  ( w  )  dw 


with 


2  exp ( A  4>  -  B  **} 

m _ m 

P(b,«j)m)  -  y  ^  eTpTT^l  -  b 


and 

Q  ( w )  -  exp  {-(A  -  2B<>  )w2  -  Bwk}. 

m 

One  can  now  optimize  P/Qdw  for  given  r,  and  T, 
its  derivative  with  respect  to  b  equal  to  zero 
gives  (after  some  algebra) 

(x  +  A ) ( x  +  a ) ( x  -  e)  +  (A  -  a)/4B  -  0, 

where 

x  -  tn(W>- 

A  -  $pk  ♦  de,/db. 


a  -  4»pk  +  e^b  +  aTt/b  -  1, 

and 

e  -  ( 2C  “ 1  -  1 ) / 4 B . 


When  is  near  Tpk  as  expected  physically,  x  is 
the  optimal  value  of  x  is 


[111. 33] 

[111. 34] 

[111. 35] 

by  setting 
;  doing  so 

[111. 36] 

[111. 37] 

[111 .38] 

[111. 39] 

[111 . 40] 
small  and 


w  aa 


[III .41  ] 


opt  4  BaA  -  N ( a  +  A ) 


with 


N  =  2C-1  -  1 . 


[Ill . 42] 


This  gives  the  optimal  value  of  b=£I2/mc2: 


„  aTpx  8  ~  aT> 

’opt  “  (e,/b)  +  (j>pk  +  x 


[111.433 


with  x  -  x  . .  Expanding  this  for  small  x,  one  has 

opt 


6opt  •  UV/A>  [’  *  ^  WrHrl-AT1!  CIII-‘,,] 


NaA  -  1 


where  the  1  dominates  the  bracket.  In  physical  terms,  if 
one  drops  the  small  correction, 


t  2  M  /  ®  \  2  (  1  2  ) 

opt  t  mi 


c2  T 


pk 


[111.45] 


or 


I  .(MA)  -  0.3  x  10*[?  (J_)t n.  (keV)A-1  ]1/2  [HI. 46] 
opt  1  £,  cm  pk 


with 


A  -  6l/b  *  ln(Tp|</Tl) 


[111.47] 


and  e,/b  is  assumed  independent  of  b.  The  radiation  inte¬ 
gral  with  the  optimal  driving  current  can  be  evaluated  using 


/ 


exp  [  -  (  C  ±  1  -2Bx  )  w  2 -Bw  -  ]dw  =  j  K  (  4,  )  (  2  B  )  ~  1  7  U 

d  1  /  *4 


4)  =  ( j  C_l  -  Bx  )  2  /  2  B 


[III. 48] 


-1/4 

which  is  approximately  0.913  when  C_1  -  2Bx  <  1/2.  This 

gives 


Y  .  ( J  )  -  0 . 2p2iA-2  l-*pj 

Anh  m  \  T  ' 


[a ( 1 +Z )C/bJ 


■g_(W  cm3) 


[r ! ( cm) /T j ( eV )  ] 


[III.  H9] 


with  u  the  total  mass  per  unit  length  (gm/cm)  and  A^  the 

atomic  mass  in  AMU.  We  note  here  that  a*A-1-3  and  A  /(1+Z) 

m 

»  2.  The  initial  radius  at  assembly,  rlf  is  expected  to  be 
of  order  0.3  cm  in  typical  experiments.  The  index  N  is 
roughly  1  / 2  . 

Before  refining  these  estimates  for  the  adiabatic  limit, 
let  us  consider  the  magnitudes  they  imply.  To  do  this,  we 
recall  section  A  of  this  chapter,  and  model  the  t<t,  cool 
collapse  phase  as  having  (1+Z)T  approximately  proportional 
to  I2(t)  once  a  Bennett  quas i -equ il i br i um  i3  reached.  The 
plasma  initially  heats  and  expands  until  reaching  such  an 
equilibrium,  although  the  equilibrium  may  not  always  be 
reached  before  time  t,.  This  collapse-phase  model  is  espe¬ 
cially  appropriate  for  wire  array  implosions  where  a  clear 
Bennett-p inch  quas 1-equlibrium  exists  for  each  wire.  Refer¬ 
ring  to  this  equilibrium  temperature  as  T0(I),  its  near¬ 
constancy  for  constant  I  gives  us  T^T,,  in  the  constant- 
current  case.  Physically,  the  near-absence  of  heating  in 
this  phase  reflects  the  fact  that  although  the  plasma  is 
accelerated  inward  by  the  field,  there  is  nothing  except 
magnetic  field  for  it  to  push  against  to  compress  and  heat, 
as  long  as  the  acceleration  timescales  are  long  compared 
with  the  sound  transit  time  across  the  plasma(s)  (i.e.,  as 
long  as  there  is  no  shock  heating). 


For  wire  array  implosions,  to  the  extent  that  the  cur¬ 
rent  I,.,  in  the  pinch  phase  is  related  to  the  typical  cur- 
rent  I  in  the  WIRES  phase,  and  to  the  extent  that  the  annu- 
lar  BRIDGE  phase  does  not  strongly  heat  or  cool  the  plasma, 
one  can  relate  r,  and  T,  at  the  beginning  of  the  pinch  phase 
to  the  driving  terms  in  the  WIRES  phase.  Over  the  limited 

range  of  total  wire  currents  I  for  which  pressur e-balance 

w 

Bennett  equilibrium  is  achieved  in  the  WIRES  phase  one  has 
for  aluminum  6-wire  arrays 

[l-t-Zdj  )  ]T1(keV)*2.25  x  1  0  “  5 1  */ u  (  g/cm  )  [III. 50] 

with  I  in  MA,  from  Equation  [III. 2].  Over  the  even-more- 

W 

limited  range  of  I  for  which  blackbody  cooling  balances 
Ohmic  heating  one  has 

r  ,  ^ ( cm  )  - 1 .5  x  10s  I ' 3Z 1 ' 3 [ {1 +Z ( T , ) } v ] 1 1 7 6  [III. 51] 

for  a  6-wire  aluminum  array  merging  at  radius  r.  .  (Again 

w 

Iw  is  in  MA,  and  the  formula  probably  does  not  apply  much 
above  I’/p-IOM  The  BRIDGE  phase  further  compresses  this 
outer  radius  by  a  factor  of  order  0.6  before  it  assembles 
into  the  pinch  phase. 

The  ratio  of  maximum  to  initial  temperatures  in  the 
pinch  phase,  T  /T,,  based  on  the  quas i -ad iaba t i c  approxima- 

l» 

tion,  is  shown  versus  the  normalized  quantity 


2 aT  , 


[III. 52] 


in  Figure  1 1 1  .  3  »  for  various  values  of  f=b/aT,.  For  alumi- 


(  I^a/u)/T, (keV) . 


[Ill .53] 
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For  any  given  generator  voltage  pulse  shape  one  expects  both 
rj  and  T,  to  be  roughly  proportional  to  b  or  Il/u,  so  that 
neither  K0  nor  r  should  vary  greatly  as  l2/p  is  changed. 
What  variation  there  is  should  come  about  because  of  changes 
in  the  wire  merge  conditions  and  duration  of  the  annular 
phase.  K0  is  ratio  of  kinetic  to  thermal  energy  content  at 
the  beginning  of  the  pinch  phase,  while  r  is  a  ratio  of 
magnetic  to  thermal  energy  at  that  time. 


2.  Singular-Perturbation  Expansion  for  Nonadiabatic  Trajec¬ 
tories 


The  energetic  radiation  output  (Equation  III. 20)  for  the 

nonadiabatic  problem  can  be  calculated  by  expanding  g(T)  and 

r(t)  about  the  peak  compression,  r-r  ,  which  we  take  for 

m 

convenience  to  occur  at  t-0.  We  keep  in  mind  the  useful 
case  where  the  nonadiabatic  trajectory  given  by 


r  -  (aT  -  b)/r 
T  -  T( -Cr/r  -  f )  , 


[III. 54] 


is  near  the  adiabatic  one  given  by  the  special  case  f-0;  in 
the  (  inr  ,  inT)  plane  the  trace  of  the  adiabatic  trajectory 
is  up  and  back  along  the  line  segment  in ( T/T , ) --C in ( r /r , )  . 
The  nonadiabatic  one  is  a  loop,  as  shown  in  Figure  III. 4, 
except  at  late  times  (when  it  may  re-compress). 

Equations  [III. 54]  give  first  integrals. 


T/T  t 


( r/r  ,  ) 


-C 


exp 


f { r ( t ) , t } dt 


[III. 55] 


and 


Cr 


-Cr  *-2a(T-T , )  +  2bln(T/T  , 


(aT-b  )  f dt . 


[III. 56] 


4 
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These  can  be  written  as 


T/T,  «■  (r/p,)"<'(1-e,(T,  T,,  a,  b,  C,  m)} 


and 


r-  ( F  -  e , ) 


1  /2 


where 


and 


..  r 

j  ti 


f  dt 


...  e -r 

J  t, 


(aT  -  b)f  dt 


and  F  is  the  function  r2  for  the  adiabatic  case. 

F  -  r2  -  (2a/C)(T-T,)  ♦  (2b/c)  ln(T/T,) 


[III. 57] 


[III. 58] 


[III .59] 


[III. 60] 


[III. 61 ] 


The  radical  in  [III. 58]  cannot  be  expanded  directly  for 
small  e2  because  F-0  at  peak  compression  in  the  adiabatic 
case.  We  will  return  to  these  formulas  when  computing  the 
values  of  r  and  T  at  peak  compression. 

For  the  radiative  yield  we  expand  the  numerator  and 
denominator  of 

f  g{x( t ) }dt 
J  r *Tt )  ; 

about  the  time  t-0  when  r  has  its  minimum  value  r  : 

m 

p(t)  ■  p  +  ?  r.t1  ♦  r  r  t1  [III. 62] 

ra  2  m  ora 

glx(t)}  -  gm  ♦  g’x  t  ♦  j(g"x2  ♦x_g')t*  ♦...[111.63] 
.  m  m  m  2mm  mm 

where  we  have  used  x  =  ln{T/T  The  subscript  m  indicates 

evaluation  of  a  quantity  at  time  0  when  r“r[n*  ( In  the  non- 
adiabatic  case  T  has  largely  reached  its  maximum  and  is 


decreasing  by  the  time  r  reaches  its  minimum  value  r  .  ) 
Figure  III. 5  shows  the  shape  1/r2(t)  for  the  adiabatic  and 
nonadiabatic  (dashed)  cases.  The  t*  term  in  r(t)  is  due 
entirely  to  the  radiative  energy  loss. 

Using  the  nonadiabatic  Equations  [111.54]  we  can  evalu¬ 
ate  the  time  derivatives  at  t-0  (r-0): 


r  »  (aT  -b)/r 
m  mm 


[III. 64] 


'r  -  aT  /r  --af(r  ,T  )T  /r  *-aC(x  )  /  ( ir  JLr  * )  [III. 65] 
m  mm  mmmm  m  m 

x  T  /T  -  -f  ( r  ,  T  )  -  -G  ( x  )/(irflr*T  )  [III. 66] 
m  mm  mm  m  mm 


x  -  T  /T  -  (T  /T  )*  ■  -C(aT  -b)r_i 
m  m  m  mm  mm 

"O'  (x  )x  (trlr*T  )“l  -  x2 
mm  mm  m 


[111.67] 


[since  T 


•C(r  /r  )T  *•  tt  (fT)|  ],  where  we  have  used  the 

rn  m  in  □  r  *  in 


fact  that  T  should  be  hot  enough  that  most  of  the  radiation 
m 

energy  loss  at  peak  compression  is  through  the  transparent 
hot  component  described  by  power  g  (  x  )  /  ( irtr  2  )  : 


T  f (r  ,T  )  -  G(x  )/US,r2) 
m  m  m  m  m 


with  G ( x  )  -  (ra/m.  )g(x)/Za 
m  i  m 


[III  .68] 


Thus  the  term  G'x  t  in  Equation  [III. 63]  i3  of  order  GG*  - 
mm 

G*x2,  the  G"x2  term  is  of  order  G*,  and  the  x  G  is  of  order 
mmmm  mm 

G*  ~  Gmx. 
m  m 

With  these  expansions  of  g(t)  and  r(t),  the  integral 
/(g/r2)dt  has  the  form 


aTm1-CBX  '*  2  {m/m,  )Z.T 

m  i  *  r 


[III  .76] 


to  first  order  in  G. 

Since  m/mi  is  the  total  number  of  ions  and  Z2  -3/2(1 +Z), 
the  ratio 


( m/m .  )  Z -T 


[III. 77] 


is  just  the  radiative  energy  loss  during  the  pinch  phase  as 
a  fraction  of  the  peak  thermal  energy  of  the  load.  Then 
Equations  [III. 75]  and  [III. 76]  can  be  used  in  [III. 73]  and 
[III. 74]  to  find 


(1  ♦  C  +  2CBX  )  v 
m  i 

(1  ♦  CBX  )  2 

m 


to  first  order  in  Y , 
[III  .78] 


i  .  e  .  i 


T_(b)  *  Tm.(b)  «xp[-(Y/2) (1 ♦C*CBX)/( 1-CBX)]  [III. 79] 


■  rmaa<6)  [III. 80] 

to  first  order  in  Y .  ^Toiad  the  maximum  T  Predicted  in 

the  adiabatic  case,  and  r  .  is  the  minimum  radius  in  the 

mad 

adiabatic  case.)  To  first  order  in  Y  then,  the  deviation 
from  the  adiabatic  trajectory  due  to  radiative  energy  loss 
reduces  the  temperature  at  maximum  compression,  for  any 
given  b«Il/ijC2,  but  does  not  appreciably  change  the  minimum 
radius  of  the  pinch.  (The  radius  change  is  second  order  in 


Subject  to  the  limitations  of  the  uniform-plasma  model 
(etc.),  the  RUNIN  code  allows  assessment  of  these  results, 
whereas  the  analytic  scaling  laws  for  the  quas i-ad iaba t i c 
radiator  do  not  apply  when  the  radiative  loss  proceeds  as 
fast  as  thermal  gain.  Because  of  this,  the  RUNIN  code  will 
shortly  be  used  for  investiga-ting  the  branching  ratio, 
Pr  ad  /  ^ d  ^  nT  ^ 1  *  as  a  funGti°n  of  *  T,  and  the  total  input 
power  -d/dt  / 1/2pV*d3r,  and  for  scoping  optimal  conditions 
at  high  currents.  Related  analytic  work  will  be  carried  on 
in  that  conjunction,  to  obtain  approximate  scaling  laws  for 
the  ef f ic  ient-rad  iator  case  as  well. 

D.  CAVEATS 

Results  such  as  these  derived  from  simplified  approxi¬ 
mate  models,  deserve  caution  in  their  application  to  real 
laboratory  systems.  We  list  here  the  major  assumptions  in 
the  foregoing  scaling  laws  and  cautions  appropriate  to 
them . 

1 .  Uniform  Plasma 

If  only  a  fraction  of  the  total  plasma  volume  radiates 
at  the  indicated  temperature,  not  only  must  the  energetic 
radiation  estimates  be  reduced  proportionally,  but  the  non- 
adiabatic  corrections  must  be  generalized  due  to  the  addi¬ 
tional  lower-temperatur e  radiation  from  the  cooler  plasma 


2.  Inefficient  Coupling  of  Kinetic  Energy  to  Radiation 

First-order  corrections  to  a  quasi-adiabatic  pinch  are 
appropriate  only  if  the  energy  losses  due  to  radiation  are 
not  extremely  large.  The  total  low-energy  radiative  loss 
from  a  cool  ohmically  heated  discharge  can  exceed  the  ther¬ 
mal  energy  content  of  the  plasma  at  any  one  time,  and  in 
fact  that  is  an  efficient  way  to  design  long-pulse  low-ener¬ 
gy-radiation  sources;  but  it  is  clearly  not  describable  by 
the  formalism  used  here  to  estimate  radiation  from  a  some¬ 
what  damped  pinch  bounce. 

3.  Energetic  Radiation  Dominates  Losses  in  the  Pinch  Phase 

•  •  ) 

The  use  of  the  approximation  Y  «  Y  in  formulas  is  not  a 
necessity,  but  is  done  for  simplicity  and  is  appropriate 
for  a  uniformly  hot  plasma  near  the  peak  of  the  radiative 
loss  function.  Formulas  [III. 78]  must  be  interpreted 
cautiously  when  applying  to  typical  laboratory  cases  where 
the  energetic  radiation  during  the  pinch  phase  may  be  a 
smaller  fraction  of  the  total  radiative  loss. 

Constant  Current 

Although  one  can  model  non-constant  currents  through 
suitable  choice  of  driver  functions  b(t)  (replacing  the 
constant  b  in  the  analysis),  we  have  not  done  so  in  this 
presentation.  As  one  can  see  from  RUNIN  code  runs,  the 
assumption  I»const.  is  only  very  approximately  satisfied 
during  the  pinch  phase. 


It' 


5.  Dependence  of  Pinch  Results  on  Initial  Conditions 

The  wires  and  annular  phases  of  the  implosion  determine 
the  initial  conditions  for  the  pinch  phase,  and  inasmuch  as 
they  depend  on  the  generator  current  (assuming  fixed  current 
versus  time  profiles  with  variable  amplitude),  the  'op t im i za- 
tions  for  pinch  radiation  done  at  fixed  initial  state  are 
not  applicable  for  the  total  problem. 

6.  Ohmic  Heating 

While  ohmic  heating  is  negligible  for  strong  implosions 

of  low-mass  load3,  high-mass  loads  accelerate  more  slowly 

and  may  have  much  lower  PdV  heating  rates.  The  smallness  of 

the  ohmic  heating  is  also  dependent  on  the  assumption  of 

fairly  complete  current  penetration  during  the  wires  or 

annular  phases ;  if  the  current  is  confined  to  a  very  thin 

sheath,  ohmic  heating  power  increases:  -  2irJ,/nJ*r  dr  ■ 

2  it  8, 1  *  /  oA  with  conductivity  o(T)  and  current-carry  ing 

curr 

area  A 

curr 

These  cautions  are  pointed  out  as  reservations  which 
could  be  major  in  some  implosions  and  some  parameter  ranges, 
but  which  are  probably  not  so  major  in  typical  present-gen¬ 
eration  experiments.  An  attempt  has  been  made  to  make  sim¬ 
plifications  to  provide  transparent  calculations,  without 
throwing  out  the  more  important  physics  affecting  the  scal¬ 
ing.  We  have  also  provided  a  quantitative  evaluation  of  one 
important  deviation  from  the  idealized  formulas  (namely  that 
due  to  radiative  cooling).  Other  deviations,  particularly 
at  the  higher  load  masses  implied  for  very-high-current 


experiments  are  less  easy  to  quantify;  for  example  some  of 
the  opacity  and  radiation  blanketing  effects  do  not  appear 
to  be  representable  by  simple  algebraic  formulas,  thus  ne¬ 
cessitating  (so  far)  case-by-case  code  calculations. 

E.  COMPARISON  OF  REQUIREMENTS  FOR  GENERATORS  DRIVING  RADIA¬ 
TION  SOURCES  AND  FUSION  TARGETS 

Inelastic  processes  leading  to  nuclear  fusion  are 
strongly  exothermic,  whereas  those  leading  to  radiation  are 
endothermic.  Because  of  this  difference,  the  nonadiabatic 
trajectories  in  lnr-lnT  space  for  radiation  and  fusion  tar¬ 
gets  curve  away  oppositely  from  the  adiabatic  trajectories. 
More  important  a  fusion  spark,  i.e.,  a  central  region  with 
enhanced  temperature,  tends  to  expand  and  propagate  the 
fusion  burn  whereas  central  radiating  spark  (enhanced  densi¬ 
ty,  or  temperature  within  limits)  tends  to  radiatively  cool 
and  allow  further  compression,  leading  to  a  radiative  pulse 
of  increasing  intensity,  but  decreasing  energy  per  photon, 
and  a  non-propagating  radiative  burn  because  of  the  high 
photon  transparency  at  low  masses.  (Thus  quasi-adiabat  ic 
calculations  tend  to  be  overly  pessimistic  for  fusion 
sources  but  overly  optimistic  for  K-line  radiation  sources. 
The  applicability  of  uniform-plasma  approximations  is  poor 
for  fusion  burn,  better  but  still  questionable  for  radia¬ 
tion  product  ion .  ) 

Another  very  important  difference  is  that  the  critical 
temperature  required  for  fusion  reactions  to  proceed  appre¬ 
ciably  is  much  higher  than  that  for  the  production  of  ener¬ 
getic  radiation  in  a  plasma  of  intermediate  z,  e.g.,  argon 
or  aluminum;  and  the  critical  fusion  temperature  is  an  ion 


temperature  whereas  the  radiation  requires  only  sufficiently 
hot  electrons  which  are  generally  easier  to  obtain.  Sub¬ 
tracting  D-T  Bremsstrahlung  energy  loss  from  D-T  fusion 
energy  gain,  the  g ( T ) -equ i valent  energy  production  curve  for 
a  fusion  neutron  source  has  a  broad  maximum  at  20-30  keV  ion 
temperature,  while  aluminum  K-line  radiation  peaks  at  1  keV 
and  argon  at  2.4  keV.  For  present  or  near-term  energy 
sources  available  to  drive  laboratory-scale  implosions,  this 
makes  radiation  sources  much  easier  to  achieve  than  fusion. 

Hydrodynamic  timescales  in  a  fusion  pellet  are  also 
faster  than  in  a  radiation  source  plasma  because  of  the 
lower  atomic  weight  of  the  fuel.  This  may  result  in  power 
delivery  requirements  for  fusion  targets  which  are  higher 
not  only  because  of  the  need  to  achieve  much  higher  tempera¬ 
tures  but  also  because  of  the  need  to  deliver  the  energy 
more  abruptly  to  avoid  pre-heating  and/or  ablating  the  fuel 
during  the  run-in  phase  of  the  liner  or  pellet. 

In  both  fusion  and  radiation  targets,  pre-heat  of  the 
central  region  during  the  run-in  phase  is  to  be  avoided 
because  it  reduces  the  peak  achievable  n2;  i.e.,  the  greater 
the  central  pressure  just  before  collapse  of  the  liner/arnu- 
lus/pellet,  the  softer  the  compression.  Since  blackbody 
radiation  from  the  imploding  plasma  annulus  can  be  a  major 
contributor  of  central  plasma  preheat  it  is  a  truism  that 
the  annulus  or  shell  should  be  kept  as  cool  as  possible 
during  the  run-in  phase.  But  in  fact  one  has  rather  little 
control  of  the  ohmic  heating  of  the  shell  or  annulus  except 
to  deliver  most  of  the  driving  power  near  the  end  of  the 
run-in.  In  laser-driven  implosions  this  is  achieved  by 
pulse  tailoring;  in  electro-magnetically  driven  implosions 


it  is  more  difficult  because  (1)  the  effective  driver  coupl¬ 
ing  is  reduced  once  the  implosion  velocity  is  high  (the 
effective  electric  field  driving  J  for  the  JxB  force  is 
E+VxB/c,  and  although  E  increases,  V  and  B  tend  to  increase 
faster  and  cancel  the  E)  and  (2)  voltage  pulse  tailoring  on 
10  ns  timescales  is  technologically  more  difficult. 


CHAPTER  IV 


SIMPLIFIED  NON-LINEAR  FIELD  DIFFUSION 


Electric  or  magnetic  field  diffusion  into  an  initially 

unmagnetized  plasma  slab  was  modeled  for  a  fixed  plasma  with 

■3/2 

temperature-dependent  conductivity  (e.g.,  *  T  )  and  with 

ohmic  heating  balanced  by  radiative  heat  loss  of  the  form 
B  (x 

RnpT  (  R ,  ct ,  B  constants).  Except  if  a  -  3/2,  the  temperature 

B  a 

is  then  given  by  JE  -  Rn  T  .  The  diffusion  equation  for  the 
electric  field  E  then  takes  the  form 

E"  -  n"pEqE  -  0, 

where  p  and  q  are  related  to  a  and  B.  This  equation  was 
solved  for  model  n(r)  profiles. 

A.  PROBLEM  DEFINITION  AND  APPROXIMATIONS 

The  diffusion  of  magnetic  field  into  a  plasma  is  affect¬ 
ed  if  the  accompanying  electric  fields  alter  the  plasma 
conductivity.  Since  ohmic  heating  raises  the  plasma  temper¬ 
ature,  the  self-consistent  field  profile  and  the  field  pene¬ 
tration  history  can  be  quite  different  from  the  ususal  dif¬ 
fusion  model,  but  can  be  calculated  if  a  suitable  equation 

for  T  is  available, 
e 

In  general,  two  effects  complicate  the  physics  of  the 
field  penetration:  fluid  motion,  i.e.,  momentum  exchange 

with  the  field;  and  thermal  conductivity,  which  makes  the 
temperature  equation  a  partial  differential  equation  as  well 
as  the  field  equation. 


But  in  certain  limiting  cases,  including  large  ion  mass, 
the  plasma  motion,  including  v  x  B  back-EMF  effects,  can  be 
ignored.  And  if  thermal  conductivity  is  sufficiently  small 
so  that  inelastic  excitation  collisions  dominate  the  heat 
loss,  we  get  a  local  ordinary  differential  equation  or  alge¬ 
braic  equation  for  T  ,  coupled  to  the  partial  differential 
equation  for  the  field. 

In  a  very  simple  case  with  fixed  ions,  qua3ineutral 
plasma,  and  negligible  thermal  conductivity,  one  may  assume 
a  quasistatic  heat  balance  between  ohmic  heating  and  excita¬ 
tion  cooling  (the  excitation  energy  is  assumed  radiated  away 
in  optically  thin  radiation): 

I  ft  (nlTl  *  Ve>  '  J'E  -  "e'xW  "o> 

where  J  -  a(n  ,  T  ) E  is  the  induced  current  density  with 
e  e 

conductivity  a  (not  necessarily  scalar),  ex  is  the  excita¬ 
tion  energy  per  ion  (or  atom),  vx,  the  excitation  rate, 
which  depends  on  the  number  density  nQ  of  ions  or  atoms  than 
can  be  excited.  If  ionization  i3  allowed  as  well  as  excita¬ 
tion,  one  has 

t  ^nt^t*ne^e*”eTe^"'r'E-neexvx^Te'no)-ne£iwi^Te,no^  CIV-2:I 

If  the  diffusion  timescale  is  slow  compared  with  the 

ionization  and  excitation  times  one  can  think  of  n^  as  de- 

e 

termined  by  T  ,  i.e.,  when 


A  V 


3  In  n 

Tt  <“.V  ■  V.('  *  i-urr  > 


CIV. 33 


with  the  logarithmic  derivative  term  a  function  of  T  .  We 


take  T 


i 


T. 


For  o(n  ,  T  )  we  have  a 
e  e 


(e*/m)  (n  /v  )  F(w  /v  ),  with 
e  m  cm 


m 


neKeTe"'3/2  +  nnKAT 
e  c  e  n  n  e 


[IV. 4] 


ui  -  eB/mc,  K'  and  K'  constants,  and  F  approaching  unity  for 
c  c  n 


The  last  term  of 
present.  For  the  fully  ionized  plasma,  a  has  the  form 


small  u)  /v 
c  m 


v  applies  if  neutrals  are 


0  *  K„TI/2  F ( <u  /neKcTe~3/2  ^  ’ 
c  e  c  e  c  e 


[IV. 5] 


on  the  other  hand,  if  neutral  collisions  dominate  \>  ,  then  a 

-1/2  m 

has  the  form  c  -  K  n  T  .  To  allow  either  case  we  let 

y  r  0  e  « 

o-Kn;T  ,  and  assume  n  <*T  n.(x)  with  some  power  e,  but 

i  e  e  e  i 


restrict  the  magnetic  field  to  values  small  enough  that 

oj  /v  <1/4  so  that  F  -  1  and  the  conductivity  i3  ap- 

c  m  - 

proximately  scalar. 


6-0 


With  the  approximation  that  oE2  and  cooling  Kxn^Te  are 


nearly  in  balance,  we  can  treat  situations  where  n  T  <<  oE: 

e  e 

and  Kxn®T“.  Then,  except  if  a  •  5,  we  have 


1 


Te  ( x  ,  t )  -  n^E2)®-5 


and 


j  ,  n^E^1 

Kx  e 


[IV. 6] 
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where  q  -  - -  and  -p  -  Y  *  5  - r. 

a-5  o-5 

(A  similar  result  {with  different  powers  p  and  q}  was  also 
obtained  for  the  "anomalously  resistive"  case  where  J/nge  is 
limited  by  the  sound  speed).  So  when  ng/ni  varies  roughly 
as  Te  we  have 


oE  +  oE  ■  Krri>  K  n  "E^E. 


This  leads  to  a  PDE  of  the  form 


EqE  -  C  nP(x)  3  ^E, 


LIV. 7] 


[IV. 8] 


(the  dot  indicates  3/3  t) ,  to  which  we  apply  boundary  condi¬ 
tions  E(-«)  *  Eq  and  E(+®)  -  0,  and  initial  conditions 


(  Eq  for  x  <  x,  <  0 
E  (x)  -  < 

*  E  exp-(x  -  x ,  )  /  A  ;  for  x  >  x. 
o  i  i 


[IV. 93 


choosing  xt  close  enough  to  where  n(x)  -  0  so  that  the  E 
wave  term  is  still  negligible  compared  with  the  oE  and/or  Eo 
terms.  A  ID  slab  model  was  used  and  n^x)  was  given  time- 
dependent  function  of  x.  in  the  wave-diffusion  equation,  Vx 
Vx  E  +  E  /c  -  -(4ir/ca)  (oE  +  Eo),  the  E  wave  term  was 
dropped.  We  assumed  one  of  two  forms  for  n^x);  a  step 
function  n,  -  niQH(x),  or  a  Guassian.,  n^^  -  n,_  exp(-x*/h*). 


n exp ( -x  */h  *  ) . 


B.  NUMERICAL  SOLUTION 

The  resulting  nonlinear  diffusion  equation  for  these 
simple  approximate  cases  was  solved  using  a  ID  finite-dif¬ 
ference  computer  code,  fully  flux-conserving,  employing  a 
fully  implicit  method  with  tridiagonal  matrix  inversion. 


The  code  was  run  for 


time  equal  to  four  initial  magnetic 


diffusion  times  for  the  Gaussian  plasma,  on  a  fifty-point 
uniform  spatial  grid.  Because  the  wave  term  was  not  kept, 


>  V'.'-ivy- ;•  ,%y  a 


WWW 


the  solutions  are  not  quantitatively  correct  far  from  the 
plasma.  The  Courant  stability  condition,  At  <  Axa/D,  was 
not  satisfied  at  the  foot  of  the  pulse  when  D(E)  had  a  pole 
at  E  -  0,  but  the  model  equations  are  not  strictly  valid  in 
that  limit  because  the  conductivity  and  temperature  do  not 
actually  go  to  zero  when  E  goes  to  zero  at  the  diffusion 
front.  On  the  timescales  of  interest,  this  failure  to  sat¬ 
isfy  the  Courant  condition  at  the  front  did  not  make  impor¬ 
tant  changes  in  the  overall  field  diffusion  history. 

C.  INTERPRETATION  OF  SOLUTIONS 

The  solutions  for  an  ionizing  Argon  plasma  with  Coulomb 

collisions  and  excitation  of  ion  lines  (B-2,  Y-0,  5-3/2)  are 

compared  with  the  solutions  of  the  standard  textbook  case,  o 

-  0,  in  Figure  IV. 1  and  IV. 2.  Where  nt  Is  not  small,  the 

buildup  of  plasma  conductivity  (for  a  >  3/2)  tends  to  short 

out  the  field  penetration  and  steepen  the  gradient  of  E.  In 

the  absence  of  other  effects,  we  would  interpret  thi3  as 

evidence  that  (1)  the  bulk  field  diffusion  in  such  a  plasma 

is  slower  than  the  rate  indicated  by  the  magnetic  diffusion 

time  4 iro  ha/ca  but  (2)  the  diffusion  of  the  "foot"  of  the 
o 

field  profile  i3  much  faster  than  for  the  constant-o  case 
because  the  conductivity  and  temperature  are  low  near  the 
foot . 

For  a  less  than  about  3  in  the  radiation/excitation  Ta 
dependence,  the  nonlinearity  of  the  diffusion,  Eq,  with  q  - 
3/(a-3/2),  can  be  quite  strong  and  sensitive  to  a;  for  Argon 
at  n  “10  cm  J  and  T  -  50-150  eV,  q  was  0.91  and  the 

I  v 

nonlinearity  was  moderate. 


\n(x) 

\ 

\ 


\  P  =  l.to 
\  Q  =  0.91 


l\2\3 


This  field  model  calculation  was  also  applied  to  a  weak¬ 
ly  ionized  region  in  a  uniform  neutral  medium,  with  elec- 

1  /  2 

tron-neutral  collisions  dominating  (Y-1,  6-  -1/2,  <ov>«T  ) 
and  with  excitation  of  neutrals  (8  -  1).  Here  p  -  -1  and  q 

-  -1/(a  +  1/2).  Again  the  diffusion  was  compared  with  the  o 

-  0  case.  Here  increasing  Tg  reduces  the  conductivity  and 
enhances  the  diffusion  rate  instead  of  impeding  diffusion, 
but  the  nonlinearity  is  weaker. 

At  higher  temperatures,  when  a  is  not  much  larger  than 
5,  i.e.,  when  the  temperature-dependence  of  the  excitation 
cooling  is  not  much  different  from  that  of  the  conductivity, 
excitation  fails  to  clamp  the  electron  temperature  at  a 
well-defined  value,  and  the  T  term  in 

'  Ve  *  G<Te)(oE2  -  Kxn®T“)  [IV. 10] 

is  not  necessarily  small  compared  with  the  right-hand  source 
and  3ink  terms.  The  equation  in  this  case  may  fail  to  spec- 

g 

ify  a  temperature  algebraically,  especially  when  o  «  n"  as 
well,  and  it  must  be  integrated  in  time  simultaneously  with 
the  E  field  equation. 

With  the  more  advanced  MHD  code  discussed  in  Part  II  of 
this  report,  we  are  also  exploring  other  effects,  including 
those  of  self-consistent  magnetization  in  the  conductivity, 
thermal  conduction,  fluid  motion,  etc.,  which  allow  a  broad¬ 
er  range  of  penetration  phenomena  for  the  electric  field. 

D.  CAVEATS  TO  THE  FOREGOING  NONLINEAR  DIFFUSION  STUD! 

To  assess  the  importance  of  ohmic  heating  and  plasma 
blowoff,  and  to  predict  the  dynamics  of  the  assembled  pinch 
plasma  if  ohmic  heating  is  important,  it  is  necessary  to 


have  some  knowledge  of  the  field  penetration  depth  into  the 
plasma.  The  process  is  clearly  affected  by  nonlinear  fea¬ 
tures  as  just  described,  but  also  by  other  features  not 
included.  Perhaps  the  most  important  of  these  is  magnetiza¬ 
tion  of  the  conductivity,  which  raises  the  effective  resis¬ 
tivity  just  behind  the  field  penetration  front  and  hence 
favors  current  flow  on  the  front  where  the  field  is  still 
low.  This  probably  tends  to  enhance  the  effective  diffu¬ 
sion  rate  of  field. 

Second,  preliminary  results  from  the  ID  code  described 
in  Part  II  of  this  report  show  that  the  plasma  motion,  not 
fully  considered  here,  is  important  not  only  in  reducing  the 
effective  electric  field  by  v  X  B  (which  i3  easily  included 
in  this  model,  taken  in  the  plasma  moving  frame)  but  also  in 
causing  compression  (nonuniform  motion)  of  the  plasma. 
Finally,  thermoelectric  fields  are  set  up,  as  shown  by  the 
ID  code  and  these  tend  to  make  a  plasma  which  is  overdriven 
(where  E  is  large)  transmit  the  force  and  cause  generator 
action  in  its  interior,  setting  up  magnetic  fields  due  to 
the  resulting  axial  current  filament.  This  may  explain  a 
feature  commonly  observed  in  the  plasma  focus,  namely  pitt¬ 
ing  of  the  cathode  on  the  axis  (presumably  by  an  ion  current 
filament);  but  such  thermoelectr ic  effects  cannot  be  includ¬ 
ed  in  the  simple  model  explored  in  this  section. 

In  addition  to  these  caveats  there  are  the  mathematical 
limitations  of  the  model.  Temperature  is  increased  by  ohmic 
heating  but  does  not  vanish  at  the  field  front  where  ohmic 
heating  vanishes.  The  radiation,  at  interesting  tempera¬ 
tures,  often  does  not  have  a  simple  power-law  fit. 

The  model  above  may  be  useful  as  a  step  in  describing 
field  penetration  in  massive  higher-Z  structures  associated 


with  plasma  implosions  but  probably  does  not  show  the  domi¬ 
nant  features  of  field  penetration  versus  field  push  in  gas 
puff  plasmas. 
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CHAPTER  V 


•-V- 


SUMMARY  OF  PART  I 


This  half  of  the  report  describes  a  simple  "scoping" 
code  for  describing  imploding  plasma  radiators,  scaling  laws 
for  energetic  radiation  from  imploded  plasmas,  and  an  effort 
to  model  the  nonlinear  field  penetration  in  the  early  stages 
of  plasma  implosion,  as  the  plasma  heats. 

RUNIN  CODE  SUMMARY 

Runs  WIRES  (optional),  then  BRIDGE  (annular  plasma), 
then  SQUEEZE  (pinch  phase). 

A.  WIRES 

1.  Evolves  position,  size,  temperature  and  current  for  N 
identical  wire  plasmas. 

2.  Pre-BBE:  wires  fatten  at  sound  speed. 

BBBE :  blackbody  (BB)  radiative  cooling  balances  ohmic 

heating,  and  with  Bennet  equilibrium  (BE),  determines 
size  and  temperature  of  individual  wire  plasmas. 

3.  Wire  centers  move  inward  via  F  -  ma ,  F  determined  by 
current  I/N  in  each  wire  ( N  wires,  total  current  I)  and 
radius  of  centers  r  .  Current  I  comes  from  given  volt¬ 
age  V(t)  via  circuit  equation  including  generator  imped¬ 
ances. 

4 .  Merge  Condition:  When  wire  plasmas  touch. 

Merge:  Replaces  wires  by  annular  shell  of  same  mass, 

volume,  temperature,  RMS  position  and  inward  velocity. 
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BRIDGE 


Evolves  outer  and  inner  radii,  temperature  and  current 
of  a  uniform  annular  plasma  shell  driven  by  V(t)  via 
circuit  equation.  Uses  variable-timestep  ODE  solver 
DGEAR  to  advance  the  variables. 

In  force  equation,  current  is  treated  as  if  acting  at 
outer  surface. 

Inner  radius  evolved  approximately  by  moving  away  from 
outer  surface  at  a  speed  related  to  the  sound  speed  (as 
one  possible  model). 

Quas i-transparent  radiative  cooling  (limited  by  black- 
body  rate),  and  ohmic  and  compressional  heating,  not 
necessarily  in  balance. 

Keeps  log  of  integrated  radiated  energy  loss  and  ener¬ 
getic  photon  loss. 

Closure  condition:  inner  radius  <  10  ^  cm. 

Closure:  3et  inner  radius  to  zero  and  go  to  SQUEEZE, 

passing  it  the  outer  raduis,  velocity,  temperature, 
current,  etc. 

SQUEEZE 

Evolves  compression  of  uniform  plasma  cylinder  driven  by 
voltage  V ( t )  and  circuit  equation,  including  quasi¬ 
transparent  radiative  loss.  (Uses  DGEAR.) 

Self-similar  motion,  driven  by  force  at  surface. 
Temperature  equation  includes  compressional,  decelera¬ 
tion  and  ohmic  heating,  and  quas i -transparent  radiative 
cooling  (limited  by  blackbody  rate).  As  elsewhere, 
ionization  state  is  controlled  instantaneously  by  tem¬ 
perature. 


4 .  Continues  log  of  integrated  radiative  energy  loss  (Y  ) 
and  energetic  photon  loss  (Y”*). 

RADIATIVE  SCALING  LAWS  SUMMARY 

Assumptions : 

1.  Almost  all  the  interesting  radiation  occurs  during  the 
pinch  phase. 

2.  The  plasma  is  uniform  (including  isothermal). 

3.  Interesting  radiation  is  transparent  and  instantaneous. 
Power  Ppad  -  irr 2 1  n2g(T),  where  n  -  ion  number  density, 
T  -  electron  and  ion  temperature,  and  g(T)  is  a  known 
function  (e.g.,  from  the  NR L  model). 

U .  Perturb  about  adiabatic  ( nonrad iat i ve  )  compression. 

5.  Initial  pinch  conditions  depend  on  history  of  the  annu¬ 
ls  plasma  and  wire  plasma  stages. 

6.  Current  I  is  used  (in  lieu  of  voltage)  as  a  scaling 
variable,  since  the  compressional  forces  depend  on  I2. 

Fundamental  Formula: 

Y  >  -  /Qn2(t)g(T(t))  ir £r  2  ( t )  d t  - 

since  n  »  ( m/m^  )  /  ( irr  2  % )  .  Here  Y>  is  the  time-integrated 

energetic-photon  yield,  l  is  the  pinch  column  length,  r(t) 
is  the  outer  radius  of  the  pinch  cylinder,  T(t)  is  the  tem¬ 
perature,  g(t)  is  the  known  radiation  function,  n(t)  is  the 
ion  number  density,  m  is  the  total  load  mass  and  m ^  is  the 
mass  of  one  ion. 


Fundamental  Approach: 

1.  Radius  r(t)  and  temperature  T(t)  are  given  to  first 
approximation  by  adiabatic  compression.  (Perturb  about 
this.)  Compression  is  driven  by  b  =  l I2/m,  assumed 
nearly  constant  during  the  pinch  phase. 

2.  Find  conditions  maximizing  Y>  and  show  how  they  scale 
with  ill2/m  and  initial  conditions  of  the  pinch,  to  low¬ 
est  order  ( quasi-adiabatically ) . 

3.  Extend  to  next  order  in  non-ad iaba t i c i ty . 

Primary  Results: 

1.  There  is  an  optimum  value  of  ll2/m  for  any  atomic  ele¬ 
ment: 


300  V(kem'  1 


-1  1 1 / 2 


[III  .46] 


where  A-e  ,/b  +  lnd  ^/T,  )  depends  on  initial  conditions, 

and  T  is  the  maximum  of  the  given  g(T)  curve  for  the 
P  K 

element.  e,  is  the  kinetic  energy  at  the  beginning  of 
the  pinch  stage,  and  is  roughly  proportional  to  b=il2/m. 
If  HI2/m  is  chosen  to  have  its  optimal  value,  Y>  can  be 
expressed  as  a  function  of  either  y=m/i  or  I.  Doing  the 
former , 


Yopt(J) 
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where  A_  is  the  atomic  mass  number  (a.m.u.),  K  is  a 
m 

dimensionless  number  of  order  unity  (see  text)  and  N  is 
a  constant  or  order  1/2,  related  to  the  adiabatic  index. 


(If  T/T,  *  (r/r,)  ,  then  N 

If  one  chooses  to  express 
than  y,  one  finds  Y>  ..  « 


radiative  loss.  Thi3  agrees  with  observations  and 

o 

simpler  scaling  models  . 

In  the  next-order  approximation,  where  radiative  energy 
loss  is  a  non-negligible  perturbation  of  adiabatic  be¬ 
havior,  falls  below  y2  or  1 4  scaling.  The  ratio  of 

Y>  to  the  peak  thermal  energy  content  of  the  load  is  a 
natural  scaling  variable  if  the  pinch  phase  is  hot 
enough  for  energetic  photons  to  dominate  the  radiative 
loss  rate.  Let  this  ratio  be  called  Y.  Then  to  first 
order  in  Y,  the  minimum  radius  (peak  compression)  of  the 
pinch  is  unaltered,  but  the  peak  temperature  is  lowered 
by  a  factor 


-  1  .  ) 

in  terms  of  I  rather 
to  zeroth  order  in  the 


exp  - ( Y / 2  )  F 


[III. 79] 


where  F  depends  on  the  adiabatic  index  C  and  on  the  peak 
temperature.  A  higher  current  is  thus  needed  to  reach 
the  optimal  radiating  temperature  than  in  the  quasi- 
adiabatic  limit,  and  this  translates  to  a  scaling  slower 
than  y 2  or  I  **  at  large  I  or  y. 


NONLINEAR  DIFFUSION  SUMMARY 


For  a  noncompr ess  in g  plasma  in  radiation  equilibrium 
.,  where  ohmic  heating  is  balanced  by  radiative  cool- 
,  we  examined  classical  field  penetration,  affected  by 


b 


the  changing  conductivity  of  the  plasma  as  it  is  heated  by 
the  penetrating  field.  This  diffusion,  in  a  medium  where 
the  diffusion  coefficient  depends  on  the  field,  shows  quan¬ 
titative  differences  from  the  usual  constant  conductivity 
field  diffusion  problem.  The  differences  depend  of  course 
on  how  the  conductivity  changes  with  field,  and  this  in  turn 
depends  on  the  radiative  properties.  Field  diffusion  in 
Argon  with  a  Gaussian  density  profile  peaking  at  10*’  - 
102ocm-J  and  with  temperatures  50  -  150  eV  was  modeled  with 
a  ID  nonlinear  diffusion  code.  Although  conductivity  build¬ 
up  in  the  plasma  somewhat  steepened  the  diffusion  front,  the 
foot  of  the  field  profile  moved  more  rapidly  than  for  the 
constant  a  case  because  the  conductivity  and  temperature  are 
low  near  the  foot.  This  study  may  have  some  applicability 
to  the  early-time  field  penetration  in  gas-puff  plasmas,  and 
to  the  location,  thickness,  and  freezing-in  of  the  current 
layer,  although  in  the  real  situation  plasma  compression 
(not  modeled  here)  is  probably  equally  important. 
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PART  II 
CHAPTER  I 

OVERVIEW  OF  THE  ELECTRODIFFUSI VE  MHD  MODEL  FOR  THE 
IMPLODING  PLASMA  RADIATOR 


Present  efforts  In  the  modeling  of  imploding  plasma 
radiators  are  centered  on  the  study  of  electrodif fusive , 
radiatively  coupled  MHD  theory  as  the  basic  tool  needed  to 
understand  the  role  played  by  each  of  the  plasma  properties 
in  shaping  the  implosion  trajectory,  emission  profile,  and 
radiation  pulse  observed  in  the  laboratory.  The  primary 
points  of  interest  in  the  plasma  behavior  are  two-tempera¬ 
ture  effects  during  compression,  chemical  potential  and 
ionization  profiles,  thermoelectric  effects  on  axis  at  as¬ 
sembly,  and  the  modes  of  current  penetration.  The  main 
points  of  interest  in  the  interaction  between  the  plasma 
load  and  the  diode  are  the  validity  of  the  diffusion  approx¬ 
imation,  the  structure  of  the  p lasma- to-vacuum  transition, 
the  fraction  of  power  absorbed  by  the  load  relative  to  that 
reflected,  and  the  resolution  in  spacetime  of  various 
sources  of  reflected  power.  Of  equal  interest  is  the  devel¬ 
opment  appropriate  computational  techniques  and  software 
with  which,  first,  to  answer  these  questions  and,  then,  to 
admit  a  smooth  extension  onto  a  much  wider  application  do¬ 
main.  In  order  to  obtain  this  kind  of  capability  it  is 
necessary  (1)  to  select  carefully  the  minimal  set  of  fields 
to  propagate,  (ii)  to  optimize  the  time  integrations  with 
respect  to  the  number  of  derivatives  needed,  and  (iii)  to 
encompass  a  wide  dynamic  range  in  number  density  and  radius. 
All  of  the  considerations  Just  set  forth  interact  and  usual¬ 
ly  conflict  to  some  degree.  Moreover,  as  shown  in  later 
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sections,  the  acute  nonlinearities  in  the  physical  theories 
both  help  and  hinder  various  aspects  of  the  task.  The  pre¬ 
sent  means  to  all  these  ends  is  the  implosion  code  ZDIPR  (Z- 
Driven  Imploding  Plasma  Radiator);  the  remainder  of  this 
chapter  is  devoted  to  a  summary  of  its  evolved  configuration 
and  performance  over  the  past  year. 

A.  ARCHITECTURE  OF  THE  IMPLOSION  CODE 

ZDIPR  may  be  rationally  subdivided  into  three  major 
subroutine  groups  -  the  physics  package  ( FLUI DOTS/TETDOTS ) , 
the  mesh  integrator  (GRIDROOT),  and  the  subcycle  integrator 
( GEARBOX/TETGEAR ) .  The  major  portion  of  the  code  and  the 
notation  appearing  in  later  discussions  have  been  documented 
in  Reference  1,  referred  to  a3  the  "original  formulation" 
for  the  remainder  of  this  report.  The  present  discussion  is 
intended  to  be  a  topical  updating  of  specific  points,  al¬ 
though  some  recapitulation  is  included. 

(1)  The  physics  package  consists  of  a  plasma  model  and  a 
field  model.  The  plasma  model  employs  both  ion  and  electron 
temperatures,  a  radial  flow  field  and  complete,  fully  magne¬ 
tized  transport  coefficients.  Theromelectr ic  fields  and 
heat  fluxes  and  appropriate  drift  speed  limits  on  J  have 
been  included  from  the  outset,  and  flux  limits  appropriate 
to  thermal  conduction  and  possible  thermoelectric  deficien¬ 
cies  have  been  added.  The  equatlon-of-state  and  radiation 
package  is  a  faithful  mimic  of  CRE  results  for  ionization 
states  and  chemical  potential,  and  utilizes  a  four  frequency 
bin  "pr obab i 1 i ty-of-e sea pe "  radiation  transport  scheme.  The 
field  model  discussed  and  explored  here  is  the  electro- 
diffusive  option,  which  is  a  good  first  approximation  when 


detailed  balance  of  incoming  and  outgoing  waves  does  not 
fail1.  The  presently  evolved  forumlation  also  has  the  sin¬ 
gular  advantage  of  checking  its  own  validity. 

(2)  The  mesh  integrator  is  now  a  new  cont inuously- 
var iable-step  predictor /corrector  accurate  through  second 

order.  Corrector  convergence  is  achieved  through  regula- 

•  8 

falsi  or  replacement  to  high  precision  (10  ).  This  purely 

iterative  convergence  scheme  eliminates  the  need  for  a  me3h 
Jacobian  and  the  attendant  machinery  for  factoring  and  sub¬ 
stituting,  in  contrast  to  the  present  subcycle  integrator. 

(3)  The  subcycle  integrator  (TETGEAR)  for  the  electro- 

* 

diffusive  option  requires  as  input,  at  various  t  in  the 
subcycle,  the  material  derivatives  generated  by  the  physics 

package  (TETDOTS),  which  in  turn  requires  the  advanced  mesh 

» 

and  velocity  fields  generated  at  t  by  the  mesh  integrator 
(GRIDROOT).  Each  derivative  call  by  the  subcycle  integrator 
thus  requires  a  full  update  of  the  mesh  to  the  intermediate 
time  point--the  primary  motivation  for  the  continuously 
variable  step  integrator. 

This  sort  of  subcycle  derives  from  the  observation  that 
a  spacetime  p.d.e.  which  has  the  derivative  operations  re¬ 
presented  by  discrete  differences  on  one  domain  of  depen¬ 
dence  appears  as  a  coupled  set  of  o.d.e.  on  the  orthogonal 
domain.  The  particular  case  here  is  that  of  spatial  differ¬ 
ential  operators,  corresponding  to  finite  difference  opera¬ 
tors  (derived  through  conventional  or  smooth  interpolants ) , 
producing  a  set  of  coupled  equations  on  the  time  domain.  In 
all  such  cases  the  coupling  of  the  time  derivates  is  expres¬ 
sed  as  a  Jacobian  matrix,  implicitly  dependent  on  the  mesh 
used  to  discretize  the  spatial  derivative  operators.  Here 
the  time  derivatives  are  the  material  derivatives  {Te, 
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0  ,  E  }  and  the  underlying  mesh  motion  {r,  r,  r,  r}  is 

6  Z 

completely  transparent  to  the  subcycling  algorithm. 

For  the  electrodif fusive  field  model  the  subcycle  vari¬ 
able  set  advanced  by  TETGEAR  consists  of  iY  ., 

E  ,  0  ...}  a  total  of  3*J  initial/boundary  value 

c  z j  c  ej  .  .  . 

problems  with  J  the  number  of  hydrodynamic  cells.  The  pre¬ 
script  ”c"  indicates  evaluation  at  the  spatial  cell  centers. 

All  external  driving  terms  enter  as  a  boundary  condition  on 

»  . 

E.  based  on  the  circuit  equation  appropriate  to  the 
®  J 

particular  diode  being  modeled.  In  the  environment  seen  by 

TETGEAR,  the  problem  is  specified  completely  by  the  input 

i  * . 

t  Y , } ,  the  interior  material  derivatives  {  Y.},  the  coupling 

J  * ,  .  J  * 

matrix  (3Y^/3Yj),  and  the  intermediate  time  points  t  . 

The  combined  fluid  evolution  subroutine  HYDROPUSH  incor¬ 
porates  the  three  foregoing  elements  and  forms  the  nucleus 
of  ZDIPR-apart  from  startup,  diagnostic  and  graphics  mod¬ 
ules.  HYDROPUSH  is  concerned  exclusively  with  the  advance 
of  <FLUI D-STATE>  over  the  variable  major  time  3tep  which  it 
selects  as  appropriate.  The  only  required  inputs  are  the 
common  blocks  <FLUID-STATE>  and  <FIELDADVDATA>  (containing 
the  previous  time  slice  of  fluid  variables  and  field  vari¬ 
ables)  and  the  appropriate  mesh  dimension,  scale  factors  and 
physical  constants,  obtained  from  a  variety  of  sources.  The 
output  is  an  update  of  <FLUID-STATE> ,  and  an  increment  of 
the  variable  MSI  (Main  Step  Index)  by  1 .  If  a  restart  file 
is  desired,  it  is  created  for  later  use  and  assigned  a  rec¬ 
ord  number  equal  to  the  MSI.  If  any  severe  errors  occur  in 
the  advance,  dump  files  are  created  to  allow  examination  of 
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intermediate  results,  and  a  variety  of  reports  at  intermedi¬ 
ate  phases  of  the  calculation  are  available  as  general 
diagnostics . 

The  sequence  of  processing  begins  by  reading  elements  of 
the  <FLUID-STATE>  into  <MESHADVDATA>  and  the  subcycle  vector 
{iY.}.  Once  the  appropriate  data  base  is  inferred  a  full 

J 

complement  of  material  derivatives  is  computed,  by  FLUIDOTS/ 
HERTZDOTS  or  in  the  electrodif fusive  mode  by  TETDOTS.  These 
material  derivatives  are  then  used  by  STEPPER  to  select  a 
major  timestep.  For  restart  purposes,  the  environment  of 
either  GEARBOX  or  TETGEAR  is  output  at  this  point  if  re¬ 
quested.  Then  a  single  subcycle  is  done  (advancing  *  Y  , 
i+1r^  and  i  +  1r^)  and  the  new  FLUID-STATE  is  output.  The 
iterative  refinement  of  the  mesh  evolution  is  now  accom¬ 
plished  within  the  calls  to  derivative  module.  The  general 
structure  of  HYDROPUSH  is  illustrated  in  Figure  1.1;  the 
details  of  the  timestep  selection  algorithms  are  discussed 
in  the  original  formulation. 

B.  MILESTONES  AND  ROADBLOCKS 

To  date  we  can  report  rather  firm  progress  in  bringing 
two  of  the  major  subroutine  groups  into  a  state  of  develop¬ 
ment  which  is  adequate  for  the  goals  set  forth  above  and  in 
the  original  formulation.  The  calculation  of  self-consis¬ 
tent  material  derivatives  is  now  routinely  done  to  very  high 
precision.  This  capability  is  documented  below  in  Chapter 
II.  A  singular  difficulty  in  obtaining  this  result  was  the 
strong  nonlinearity  in  the  Br=ginskii  thermoelectric  theory, 
but  this  has  now  been  tamed  and  at  least  partially  under¬ 
stood  as  a  possible  basis  for  corrections  to  that  theory. 


FLOW  OF  A  SINGLE  SUBCYCLE:  HYDROPUSH 


Input:  V  and  [T  E_  0  ]  *  Y  at  I-th  time  level 
r  I  Z  e 


Set  major  timestep:  CFL,  CELL  AREA  CHANGES,  SURG 


Output  RESTART  files,  if  requested 


ENTER  SUBCYCLE  INTEGRATOR 

-[A]  Set  minor  timestep,  t^ 

PUSH  PREDICTOR:  Y . 

1+1 10 

CONVERGE  CORRECTOR:  Y.  „(Y.  ,) 

1+1  1+1 


-[B]  Select  JACOBIAN,  call  TETJAC  or  keep  old  one 

Calculate  iterated  derivatives  Y.  ,  (TETDOTS) 

i+l  ,m 


L*-Test  step,  order  &  error,  recycle  as  necessary 
Recycle  over  minor  timestep  until  major  timestep  is  reached 

Output  the  advanced  set  to  the  I+l  time  level,  update  r,  history. 


Figure  1.1 
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Also  completely  operational  is  the  mesh  integrator,  documen¬ 
ted  in  Chapter  III.  The  development  of  these  new  integra¬ 
tion  filters  grew  from  the  great  sensitivity  of  the  subcycle 
integration  to  the  space-time  jitter  in  the  velocity  field. 
Their  subsequent  introduction  in  a  fully  implicit  implemen¬ 
tation  has  improved  the  quality  of  the  model  significantly. 

The  primary  obstacle  to  efficient  operation  is  the  sub¬ 
cycle  integrator  DGEAR,  essentially  unchanged  from  the  ver¬ 
sion  supplied  by  IMSL.  While  the  integrator  DGEAR  offers 
several  control  parameters  and  considerable  flexibility  in 
the  integration  scheme,  the  most  useful  options  have  been 
the  use  of  explicit  Jacobians  and  the  control  of  the  local 
truncation  error  through  a  variable  TOL.  The  explicit  Jaco¬ 
bian  is  the  most  favorable  option  in  this  hydrodynamic  ap¬ 
plication;  but,  since  the  physics  detail  in  ZDIPR  demands  a 
very  expensive  analytical  evaluation  of  this  matrix,  an 
estimate  based  only  on  finite  differences  is  needed.  We 
find  that  the  choice  of  the  differencing  algorithm  is  a 
sensitive  one  relative  to  subcycle  integrator  performance. 
The  variable  TOL  provides  a  natural  and  consistent  means  of 
specifying  time  integration  parameters,  but  it  does  not 
split  up  into  separate  error  bounds  for  truncation  and  cor¬ 
rector,  and  this  causes  trouble.  The  integration  scheme  can 
be  selected  as  either  an  implicit  Adams  method  of  up  to 
twelfth  order  or  a  backward  differentiation  method  of  up  to 
fifth  order  (Gear’s  stiff  method).  Both  methods  are  of  the 
implicit  linear  multilevel  type  and  require  the  solution  of 
an  algebraic  system  at  each  interior  (subcycle)  timestep. 
Here  the  optimal  choice  between  the  two  seems  to  fall  to 
Gear's  stiff  method.  The  motivation  for  its  use-treatment 
of  stiff  terms  in  the  heat  and  field  diffusion  equations- 
therefore  appears  well  grounded;  but  the  error  estimates  the 


step  size  controls,  and  the  dependence  on  a  Jacobian  matrix 
(and  its  inverse)  for  corrector  convergence  appear  to  con¬ 
tain  some  systematic  and  serious  difficulties  in  this  appli¬ 
cation.  The  problems  are  being  addressed  by  a  variety  of 
alternatives  and  the  resolution  of  these  questions  will  be 
best  assured  by  careful  further  testing. 

Despite  the  difficulties  in  obtaining  efficient  integra¬ 
tion,  the  model  can  be  fruitfully  applied  in  short  simula¬ 
tions.  Chapter  IV  describes  the  results  of  one  such  calcu¬ 
lation  meant  to  examine  the  early  boost  phase  of  the  plasma 
implosion . 


CHAPTER  II 


SELF-CONSISTENT  CALCULATION  OF  THE  MATERIAL  DERIVATIVES 

( FLUIDOTS/TETDOTS) 


« 
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A  number  of  modifications  to  the  original  formulation 

g 

have  been  necessary  in  the  course  of  the  code  development  r 

r 

over  this  contract  year.  These  changes  have  been  motivated  £ 

primarily  by  the  need  for  high  precision  in  the  difference- 
equation  image  of  the  time  derivatives.  The  better  resolu-  i 

tion  of  the  global  couplings  (in  space)  for  each  temporally  ® 

evolved  quantity,  as  provided  by  this  high  precision,  has  <’ 

been  found  very  helpful  in  improving  the  subcycle 
integration. 

A.  STATEMENT  OF  THE  PROBLEM 

The  physical  content  of  this  model  has  been  laid  out  in 
the  original  formulation  and  in  previous  documents.  How¬ 
ever,  in  order  to  clarify  the  present  discussion,  the  major  ^ 

** 

fluid  and  field  evolution  relations  are  summarized  below  for 

s' 

the  electrodif  f  us  ive  limit  that  has  been  the  central  focus  jj 

of  recent  efforts.  The  notation  is  unchanged  from  that  H 

introduced  in  the  original  formulation.  > 

If  one  gives  up  some  information  concerning  the  details 
of  the  diode  fields  and  makes  the  assumption  that  the  incom-  J 

ing  and  outgoing  wave  components  are  in  detailed  balance,  s; 

then  the  Hertz  wave  equation  can  be  transformed  to  a  diffu- 
sion  equation  ;! 
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where  E.h  is  the  dimensionless  thermoelectric  field,  Z  is  a 
t  n  a 

dimensionless  conductivity,  8«Vfluid/c,  E'-E  +  BB,E  »  E  + 

BB  +  E..  ,  and  D/Dt  -3  +  83  . 
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The  fluid  response  to  the  electromagnetic  stresses  and 
heating  is  embodied  in  the  relations 
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In  these  expressions  9 
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X  T  is  the  thermal  conductivity,  t  the  plasma 
6  )  1  6 


relaxation  time,  Er  the  ambipolar  radial  field  (with  p  its 


rad 


is  the  net  (local) 


induced  charge  density),  and 
radiative  heating  or  cooling.  The  dimensional  version  of 
fields  are  subscripted  with  a  vector  component;  dimension¬ 
less  fields  are  not.  The  radial  electric  field  is  estab¬ 
lished  as  a  solution  to  an  integral  equation  derived  from 
the  radial  component  of  Ohm’s  law.  The  dr  if t-speed- 1 imi ted 
current  condition  is  supplemented  by  a  (nonlinear)  change  in 
Z  where  the  local  E  field  requires  it. 
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The  (nonuniform)  time  levels  for  any  variable  are  index¬ 
ed  with  a  leading  superscript;  the  Lagrangian  fluid  mesh  is 
denoted  by  r,  and  its  material  derivatives  by  a  superscript¬ 
ed  dot  (or  dots).  Spatial  indexing  is  denoted  by  trailing 
subscript  and  various  cell-to-cell  averaging  operations  are 
denoted  by  an  overbar  or  by  angle  brackets. 

The  data  base  for  all  the  hydrodynamic  calculations  is 

,  c 
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the  common  block  <FLUID-STATE> ,  containing  {N.,  ,Tt  .  , 

i  ii*  i  J  1  L * J 

e  .,  r  ,  r.}  and  {  E_  .},  the  diffusing  electric  field. 

ce,jjj  cz,j 

An  update  of  the  < FLUI D-STATE >  and  Ez  is  the  central  result 
of  a  major  timestep.  The  relationships  among  the  basic 
fluid  variables  on  the  Lagrangian  mesh  are  illustrated  in 
Figure  II. 1.  The  simple  two  and  three  point  area-weighted 


differencing  schemes  for  such 


mesh  were  discussed  pre¬ 


viously.  The  numbers  of  cell  ions  { N  .  }  are  a  conserved 

J  i 

vector  of  ions/cm  resident  in  the  (compressible)  cell  [  r,, 

i  ^ 

rj+1],  assuring  strict  particle  conservation  and  a  solution 

of  the  equation  of  continuity  limited  in  accuracy  only  by 


the  evolved  values  of  (  r j } . 


The  elements  of  { N ^ }  are  as¬ 


signed  spatial  locations  given  by  the  cell  center  position 
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(defined  by  the  equal  area  point) 
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In  general,  at  any  point 


(selected  by  the  subcycle 


integrator)  within  the  subcycle  ^lo’^hi^*  °ne  must  maP 
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intermediate  fields  {  r 


Vrj  cB0j 


are  all 


consistent  with  one  another.  Only  the  time  can  enter  as  a 


generator  of 


however  these  fields  are 
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calculated,  because  of  the  Lagrangian  nature  of  the  mesh. 
The  calculation  of  these  position  and  velocity  coordinates 
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considered  given  for  the  remainder  of  this  discussion. 
*  * 

Although  the  V  .  and  B„.  calculations  must  be  considered 
rj  c  6j 

coupled  problems,  a  magnetic  field  calculation  from  any 
given  velocity  profile  is  the  pivotal  element  in  achieving 
the  desired  self-consistency. 


Notation 

Figure  II. 1.  The  Lagrangian  Fluid  Mesh. 


B.  THE  MAGNETIC  FIELD  COMPUTATION 

If  the  mesh  and  its  velocity  are  presumed  given,  or  at 
least  provisionally  established  subject  to  iterative  refine¬ 
ment,  then  there  is  just  one  magnetic  field  profile  (and  its 


curl,  the  current  density)  which  is  compatible  with  this 

mesh  and  the  input  electric  and  temperature  fields.  Because 

of  the  great  sensitivity  of  the  transport  coefficients  (in 

this  often  highly  magnetized  model  plasma)  to  the  value  of 
* 

B..,  the  stability  and  accuracy  of  the. fluid  integration  i3 
c  b  j  # 

improved  considerbly  as  the  resolution  of  cB0j  is  improved. 
This  calculation  however  always  introduces  a  nonlocal 
constraint  into  the  model. 

For  the  plane  parallel  waveguide,  one  may  derive  either 
electric  or  magnetic  diffusion  relationships  from  the  Hertz 
potential  wave  equation  in  the  limit  of  detailed  balance 
between  incoming  and  outgoing  waves.  This  limit  is  also 
equivalent  to  the  neglect  of  displacement  current  relative 
to  conduction  current.  Because  both  E  and  B  are  solenoidal, 
the  resolution  of  the  "companion"  field  (E  when  B  is  diffus¬ 
ed,  B  when  E  is  diffused)  always  involves  the  calculation  of 
the  field  from  its  curl.  It  is  this  calculation  which  pro¬ 
duces  the  nonlocal  constraints  that  supplement  the  local 
diffusive  evolution  relationships.  Whether  one  chooses 
electric  or  magnetic  diffusion,  it  is  easy  to  see  that  this 
nonlocal  constraint  is  always  equivalent  to  a  single  inte¬ 
gral  equation. 


Turning  first  to  magnetic  diffusion,  the  evolution  rela¬ 


tion 
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is  obtained  by  recovering  the  curl  of  Ampere's  law  from  the 

A 

Hertz  potential  wave  equation,  writing  J-EBr,ag  insj<  i  ^  E*  and 
eliminating  3^E.  The  fields  B,  8p,  E,  radial  coordinate  u, 
time  t,  and  conductivity  I  are  in  the  dimensionless 
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form  of  the  original  formulation.  The  contribution  required 
from  Braginskii's  kinetic  theory  and  transport  coefficients 
is  the  const ituit ive  relation 


Jz  .  t(Ez.  8rBe  .  Eth>  a  £  E  ,  [II. 5) 

where  I  and  Efch  are  functions  of  BQ  and  the  plasma  relaxa¬ 
tion  time,  Tg.  On  its  face  [II.  4]  is  a  fully  local  rela¬ 
tionship,  but  one  needs  a  complete  profile  E  (u)  given  at 

z 

all  times  in  order  to  calculate  the  material  derivative  on 

the  left.  Since  B.  and  3tBb  are  the  only  fields  available 

one  must  use  Faraday's  law  to  obtain  E  (u)-/gdQ(3  B_),  a 

nonlocal  operation.  This  forces  one  to  solve  an  integral 

equation  for  D/D  B  ,  viz.  [II.  U]  is  now  an  integral  equa- 

T  0 

tion  for  the  spatial  profile  of  the  material  derivative. 
The  value  of  D/D  Bq  at  u»u,  depends  upon  the  values  estab¬ 
lished  on  the  domain  (o,  u,)  through  the  field  Ez(u).  Be¬ 

cause  [II.  i*]  is  in  fact  equivalent  to  Ampere’s  law,  however, 
to  solve  it  as  an  integral  equation  is  to  establish  that 
profile  J  which  satisfies  the  const ituitive  relation 

Z 

[II. 5].  The  only  difference  is  that  it  is  D/D  B  ,  rather 

T  0 

than  J  itself,  which  is  adjusted  to  obtain  agreement.  One 

Z 

can  equally  well  view  [II. 5]  therefore  as  the  Integral  equa¬ 
tion  being  solved  in  a  magnetic  diffusion  calculation;  the 
two  are  equivalent,  interchangeable  problems. 

For  electric  diffusion,  on  the  other  hand,  the  evolution 
relationship  is  derived  by  applying  the  Laplacian  to  the 
Hertz  potential  wave  equation,  eliminating  the  3  *72  z 
contribution,  and  introducing  [II. 5]  as  the  const itu it i ve 
relation.  The  evolution  relation  obtained, 

ffr  V  inq  (u',9u  u3ue2>-K‘"  1  •  Vth-VxV  ClI-6] 
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is  equivalent  to  [II. 4]  in  all  respects.  The  nonlocal  con¬ 
straint  however  now  enters  directly  in  the  form  [II. 5]  be¬ 
cause  the  additional  needed  profile  is  now  Ba  which  must  be 
evaluated  from  its  curl,  4irJ,  so  one  must  solve  [II. 5]  as  an 
integral  equation  in  order  to  establsh  the  profile  DEz/D  . 

Having  established  [II. 5]  as  the  needed  integral  equa¬ 
tion  for  either  form  of  the  field  diffusion  theory,  it  is 
natural  to  ask  whether  or  not  a  particular  const i tu i t i ve 
relation,  in  this  case  Brag  ins k ii ' s ,  poses  a  soluble  problem 
in  the  practical  context  of  a  known  discrete  mesh,  and  given 
discrete  images  of  the  input  fields.  If  an  accurate, 
unique,  and  continuous  solution  cannot  be  established  on  a 
given  mesh  then  the  cons t i tu i t i ve  relation  is  possibly  in 
error  -  for  what  [II. 5]  establishes  is  that  a  B„  field  and 
its  curl,  J,  are  related  in  a  particular  way.  A  result 
well-known  from  first  principles  is  that  these  fields  are 
related  uniquely  whatever  the  particular  details  imposed  by 
the  medium;  any  ambiguity  in  the  solution  raises  the  choice 
between  giving  up  on  a  Lagrangian  formulation  or  questioning 
the  accuracy  of  the  cons t i tu i t i ve  relation. 

Using  Braginskii's  notation  (cf.  Reference  2,  p.  268 
(Eq.  6.18-6.20)  and  p.  249  (Eq.  4.32  and  4.35))  and  casting 
the  remaining  discussion  in  CGS  units,  the  axial  components 
of  the  "Ohm's  law"  comprising  [II. 5]  can  be  written 

J2  ■  v* 
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with  the  further  reduction 
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obtained  when  all  the  substitutions  are  carried  out.  The 
coefficients  B’J ,  B? ,  <5  0  ,  6  t  are  functions  of  the  ionization 

A 

state  only.  When  E  ,  J_.  and  B.  are  all  negative  (the 

Z  Z  0 

usual  convention  in  ZDIPR)  or  all  positive,  and  B  is 

r 

negative,  the  domains  of  positive  T  weaken  the  net  field 

r*  © 

by  opposing  the  applied  E  .  The  domains  of  negative  T  , 

z  r  6 

in  contrast,  strengthen  the  net  field  because  Efch  is  paral¬ 
lel  to  E  .  In  this  way  the  E..  tends  to  oppose,  through  the 
z  .  tn 

J2  component  it  induces,  changes  in  the  flow  Bp  arising  from 
the  temperature  gradient. 

Using  these  results  the  basic  constraint  relation  [II. 5] 
becomes 
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[II. 9] 


with  B  -  2/cr  Jr0{2nrlJ  )  dr1P  and  x  -  -(1.7588*10T  t.)  B,, 
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making  manifest  its  characterization  as  a  nonlinear  integral 
equation . 


w 


C.  APPLICATION  OF  THE  CONSTRAINT  RELATION 


The  nonlocal  constraint  requires  in  many  cases  some  new 
modifications  to  the  field  diffusion  due  to  large  and  quite 
possibly  unphy3ical  thermoelectric  fluxes. 


In  a  spatially  inhomogeneous  raagnetoplasraa  the  process 
of  radial  implosion  forced  by  axial  currents  produces  radial 
temperature  gradients.  These  gradients  feedback  either 
positively  or  negatively  on  the  originating  axial  current 
(J_)  by  means  of  thermoelectric  momentum  fluxes  arising  in 
the  magnetized  portions  of  the  annular  load.  All  such  ther¬ 
moelectric  momentum  transfers  originate  in  differing  mean 
rates  of  electron  (r-z)  gyration  (in  the  local  magnetic 
field)  on  the  "up  hill"  and  "down  hill"  sides  of  any  phase 
space  point.  As  noted  above,  the  explicit  calculation  of 
these  imbalanced  mean  fluxes  constitutes  a  major  contribu¬ 
tion  of  plasma  kinetic  theory  to  the  expected  response  of 
the  fluid  Jz ( Ez » Bg  ,  Bp  ,Tg )  in  applied  electromagnetic  fields. 

The  constraint  [II. 9]  discussed  above  can  be  viewed  as  a 
check  or  any  constituitive  relation  proposed  from  kinetic 
theory  (or  another  source).  If  one  fails  to  find  a  continu¬ 
ous  and  single-valued  solution,  then  the  proposed  response 
Jz(Ez’B9’0r’V  is  3usPect  or  non-integrable .  In  the  course 
of  testing  the  present  ID  implosion  code  (ZDIPR)  the  inte¬ 
gral  equation  for  J  implied  by  Braginskii  theory  was  ex- 
amlned  in  this  light  and  found  to  cause  integrating  failures 
when  large  temperature  gradients  were  encountered  and  when 

thermoelectric  fields  began  to  add  to  the  applied  E  .  These 

z 

are  unfortunately  conditions  which  are  quite  common  in  the 
implosion  solutions  examined  to  date,  and  it  was  thus  neces¬ 
sary  to  find  a  remedy  if  any  progress  was  to  be  made  in  the 
physical  understanding  of  the  current  penetration  process. 
The  failure  in  integrating  Braginskii's  theory  manifests 
itself  in  a  mult  1  valued  solution  of  the  J  (r,t)  problem  when 

Z 

given  a  set  of  profiles  (o(r,t),  8r(r,t),  T  (r,t)  and 

E  (r,t)). 

Z 


This  must  be  rejected  immediately  as  a  violation  of 

uniqueness  and  a  closer  examination  shows  that  continuity  in 

the  solution  is  failing  as  well.  The  net  result  is  that  for 

a  given  non-uniform  fluid  mesh  there  always  exists  an  upper 

limit  to  the  value  of  la  T  I  that  can  be  admitted  in  calcu- 

i  re' 

lating  the  thermoelectric  fields,  Eth’  by  any  implicit 

scheme  over  the  spatial  mesh.  Somewhat  les3  obvious  is  the 

result  that  this  upper  limit  depends  on  the  radial  flow  Sp — 

Braginskii's  kinetic  theory  treatment  completely  decouples 

B_  and  E.  .  . 
r  t  h 

It  is  also  easy  to  verify  that  the  usual  validity  condi¬ 
tions  on  Braginskii’s  theory  are  not  violated  when  this  kind 


of  failure  arises. 


the  thermoelectric  field  was  limited  to 


For  example,  In  one  ZDIPR  calculation 
eld  was  limited  to  -0.03  •  E..  ,  Bra¬ 


ginskii's  on  a  weakly  magnetized  (w  t  -10  )  interior  do- 

©  6 

main  of  the  mesh  in  order  to  avoid  multivalued  roots  for  J  . 

z 

On  this  same  domain,  the  Coulomb  logarithm  A-9.87;  the  ther¬ 
mal  gradient  scale  length  L-0.2194  (cm)  is  large  compared  to 
the  mean  free  path  4-4.795  x  10  ^  (cm);  and  the  time  scale 
for  temperature  change  t-3.356  (ns)  is  large  compared  to  the 

......  M.l.u.fct.a  .  I  _  _  )■  *5  O  1  _  1  ft  3  /  \  rv.  ....  k.. 


electron  relaxation  time  t  -4.221  •  10  (ns).  On  another, 

e 

highly  magnetized,  exterior  domain,  the  thermoelectric  field 

was  limited  to  -0.0023  E..  4 .  in  the  presence  of  a 

th  Braginskii 

less  stiff  temperature  gradient  L.  -0.3518  (cm),  and  a  very 

-  4 

small  gyroradiu3  rg-.84  •  10  (cm).  In  short  it  would  not 

be  possible  to  anticipate  this  failure  on  the  basis  of  the 
constraints  on  validity  given  by  Braginskii.  Even  though 
similar  kinds  of  limits  to  7T@  apply  in  thermal  conduction, 
the  problem  here  cannot  be  formulated  with  reference  to  only 
one  process,  such  as  the  formation  of  a  heat  flux  in  the 
presence  of  a  given  magnetic  field.  Here  the  difficulty  is 
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not  local,  but  global,  and  arises  only  when  the  fluid  and 
field  models  are  coupled. 

This  limit  on  the  practical  integrability  (to  some  reli¬ 
able  accuracy)  of  Braginskii's  theory  is  a  serious  obstacle 
to  its  inclusion  in  the  implosion  model.  If  one  accepts  the 
theory  as  valid  and  resorts  to  an  explicit  integrator,  then 
the  accuracy  will  suffer  and  color  the  results  in  a  very 
subtle  way.  If  one  modifies  the  Braginskii  picture  explic¬ 
itly  to  prevent  the  problem  and  to  allow  the  more  accurate 
implicit  integrator  then  the  physical  content  of  the  model 
may  be  hurt  but  the  extent  and  degree  of  the  inaccuracy  is 
known.  The  source  of  the  difficulty  is  the  very  positive 
nonlinear  feedback  of  the  thermoelectric  field  on  itself, 
and  it  is  worth  asking  if  this  should  be  so. 

If  the  nonlinear  Braginskii  picture  of  the  thermoelec¬ 
tric  effect  is  in  error  as  the  temperature  gradient  in¬ 
creases  relative  to  B  ,  the  correct  expression  for  E..  must 

p  v  n 

contain  the  same  symmetries,  however  it  is  derived.  Since 
the  primary  source  of  this  effect  is  the  thermal  gradient, 
one  solution  is  to  limit  the  momentum  flux  produced  by  this 
gradient  while  not  altering  the  gradient  per  se,  in  common 
with  similar  treatments  for  heat  conduction.  Because  the 
magnetic  field  serves  only  to  redirect  the  thermally  gene¬ 
rated  momentum  flux  and  to  limit  the  phase  space  volume 
contributing  to  the  mean  flux  (sought  as  a  source  for  Efch), 
a  plausible  model  would  be  to  replace 
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with 


th 


3rTe)  8.(x), 


where  F  is  an  odd,  saturating  function  of  We  3^Te  (the 
basic  measure  of  source  strength)  and  x-wg rg .  In  general  F 
must  couple  to  the  flow  field  8r  and  to  the  external  field 
E  .  The  product  of  separately  odd  functions,  F(l/eVT  )  and 

Z  6 

8~(x),  is  preserved  by  the  cross  product  nature  of  the 
transport . 

The  general  expectation  one  has  for  F  is  that  it  should 

be  linear  for  small  values  of  the  argument,  and  bounded 

above  for  those  situations  where  1 /e  3  T  — >  «.  The  value 

r  e 

F ( A )  that  constitutes  this  upper  bound  can  be  set  from  the 

requirement  that  no  multivalued  roots  occur  in  i- he  local 

response  curve:  J  -  o(E„  ♦  8  Ba  +  E.. ).  In  general  one 

z  z  r  8  th 

would  expect  modifications  in  8~(x)  also,  but  the  qualita¬ 
tive  form  should  be  the  same.  At  least  the  8~(x)  dependence 
calculated  by  Braginskii  provides  a  reasonably  straightfor¬ 
ward  method  for  limiting  F(A). 

A  complete  and  necessarily  elaborate  investigation  of 
this  limiting  effect  is  certainly  needed,  but  it  is  a  very 
large  effort.  In  order  to  achieve  the  new  required  coupl¬ 
ings  (to  the  mean  flow  8  and  the  imposed  E  )  one  must  keep 

*  z 

more  terms  in  the  collision  integral,  retain  some  nonlinear 
features,  and  probably  include  density  gradients  as  well. 
Such  a  calculation  is  most  fruitful  on  a  large  scale,  ob¬ 
taining  corrections  to  all  the  transport  processes.  The 

choice  here,  to  limit  the  effective  magnitude  of  3  T  ,  is  ad 

r  e 

hoc  but  physically  and  numerically  reasonable  in  the  face  of 
such  gross  violations  of  the  electrodynamic  constraints 
discussed  above.  Moreover,  this  limit  can  be  incorporated 


smoothly  into  the  other  thermoelectric  terms,  viz.  the 
thermoelectric  heat  flux  and  the  3  E.  .  required  in  the  elec- 

x  v  n 

trie  diffusion  equation,  and  it  can  be  extended  easily  to  2D 
models . 

The  solution  of  the  integral  Equation  [II.  9]  proceeds 
from  the  mesh  interior  outward,  seeking  to  establish,  from 
the  current  densities  already  known,  a  local  current  densi¬ 
ty,  c,  as  the  next  additional  element  in  the  solution 
vector  . 

The  limit  is  derived  by  constraining  the 

3SP"dI(Ez'  V  Br’  V  <  '■ 

where  c  is  the  added  local  current  density  at  any  spatial 
point,  and  J  (E  ,  B  ,  8  ,  T  )  is  the  response  curve  given  by 

Z  Z  0  P  © 

the  transport  theory.  The  saturated  amplitude  for  F(1/eVT  ) 

is  replaced  by  (A/e)  3  rTe  and  a  value  for  A  is  calculated 

from  the  constraint  on  J'.  If  this  A  value  is  in  [0,1]  then 

z  . 

the  thermal  gradient  3rTe  is  limited  to  A  3pTe  in  that 

spatial  domain  for  all  thermoelectric  effects.  Limits  on 

the  colli3ional  thermal  flux  are  treated  separately  in  a 

similar  manner,  and  employ  a  different  limit. 

Once  implemented,  one  finds  that  this  limit  arises  under 

conditions  quite  easy  to  interpret  physically:  (1)  the 

( ExB )  drift  speed  (from  E..  xB_)  in  the  given  thermal 

p  tnfz  y 

gradient  must  tend  to  exceed  the  radial  flow  8  :  or  (2)  the 

r 

thermoelectric  ExB  drift  speed  must  exceed  the  local  dif¬ 
fusion  speed  of  B  ;  and  finally,  (3)  the  thermoelectric 

field  must  produce  a  posit  i  ve  feedback  on  the  current  in  the 
plasma.  This  is  very  much  the  situation  upon  closure  of  the 
fluid  annulus  for  example.  Particularly  interesting  among 
these  criteria  for  the  flux  limit  A  is  the  involvement  of 


V  ^  ^  / 


the  new  E  x  B  drift  speed.  It  is  therefore  at  least  a  plau¬ 
sible  conjecture  that  Braginskii's  theory  produces  this 
extreme  and  perhaps  incorrect  nonlinearity  because  the 
drift  is  not  properly  included  in  the  orbit  integrals  of  the 
Landau  collision  term  —  this  collision  term  in  fact  neglects 
all  external  fields. 

These  modifications  to  have  a  significant  impact  on 

the  physical  predictions  of  the  ZDIPR  model  as  well  as  the 
numerical  techniques  required  for  either  electric  or  magnet¬ 
ic  diffusion.  The  most  significant  physical  result  is  the 
level  of  enhanced  current  predicted  on  axis  when  the  annulus 
closes  at  high  velocity.  Under  these  conditions  the  (per¬ 
haps  small)  initial  value  of  Ba  in  the  plasma  interior  is 
amplified  by  the  thermoelectric  effect  because  the  tempera¬ 
ture  gradient  is  negative  and,  since  Ez  and  Efch  are  there¬ 
fore  parallel,  a  positive  feedback  on  J  is  generatd  by  E  . 

z  cn  . 

This  produces  no  more  J  •  E  heating  because  the  thermoelec¬ 
tric  heat  flux  will  tend  to  cool  the  axis,  but  it  does  tend 
to  increase  B.  near  the  axis  and  accrete  incoming  plasma  on 
a  radially  growing  central  filament  because  the  resulting 
J z B g  stress  rises  to  overcome  the  outward  pressure  gradient 
force,  which  might  ordinarily  slow  or  reflect  the  incoming 
plasma.  Of  course  very  large  currents  on  axis  tend  to  mag¬ 
netically  insulate  those  imploding  regions  just  outside  and 

soften  the  effect.  In  ZDIPR  test  runs  with  the  modified  E.. 

th 

the  axial  current  peak  was  still  quite  high  but  the  cells 

outside  the  minimum  in  T  did  not  tend  to  insulate  as  much 

e 

as  before  the  effective  temperature  gradient  was  limited. 
This  produced  a  more  persistent,  narrower,  and  stable  peak 
feature  in  the  current  density  profile  and  number  density 
profile,  indicating  a  (quite  reasonable)  coupled  compression 
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of  particle  density  and  current  density.  The  thermal  cur¬ 
rent  source  simply  tries  to  preserve  the  ratio  BQ/n  as  the 
plasma  stagnates,  a  result  familiar  from  magnetic  diffusion 
theory  in  the  limit  of  high  conductivity  quite  appropriate 
to  the  hot  axial  plasma.  Alternately  viewed,  the  source 
term  *  -B  p0p  in  DBQ/Dt,  tends  to  dominate  the  hot,  stagnat¬ 
ing  plasma  and  requires  current  density  to  rise  inside  any 
given  radius  in  order  to  increase  the  BQ  value  at  that  radi- 
us.  In  an  electric  diffusion  solution  this  is  done  only  by 
generating  a  strong  enough  electric  field  to  support  the 
needed  Jz>  and  the  uhermoelectric  field  (properly  limited  as 
shown  above)  is  the  source  of  this  current  until  the  corre¬ 
sponding  sources  in  the  electric  diffusion  equation,  viz. 
-  3  E..  and  -B 9  0  can  build  up  the  "external"  field  suffi- 

Toil  T 

ciently  to  do  the  job.  Such  high  axial  current  densities 
are  likely  to  produce  anode  damage  in  small  areas  corre¬ 
sponding  to  the  spatial  extent  of  the  central  current  fila¬ 
ment.  Their  evolution  is  therfore  a  high  priority  in  future 
work . 

A  second  physical  result  occurs  at  the  exterior  edge  of 
the  annulus  where  the  thermoelectric  field  tends  to  cancel, 
rather  than  to  enhance,  the  applied  E.  Here  the  limit  on 
3pT e  is  invoked  during  compression  of  the  outer  layers  with 
A  [0.800  —>  0.900],  This  is  a  consequence  of  the  required 
odd  parity  of  F  (1/e  9pTg)  discussed  above,  because  there  is 
really  no  danger  of  multiple  roots  when  a^T^O  and  0p<O. 
The  limit  here  forces  the  plasma  model  to  be  less  effective 
in  shielding  out  the  applied  Ez  and  can  be  expected  to  alter 
the  all-important  field  penetration  process  when  integrated 
over  the  run-in  phase  of  a  typical  implosion. 
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D.  PROCESSING  SEQUENCE  FOR  THE  MATERIAL  DERIVATIVES: 


FLUIDOTS/TETDOTS 

When  one  combines  the  foregoing  results  with  the  re¬ 
quirements  stemming  from  the  basic  algebraic  structure  of 
the  derivatives,  there  remains  very  little  latitude  in  the 

processing  sequence  required.  The  first  task  is  to  estab- 

*  * 

lish  the  new  mesh  {  r.,  V  . }  and  the  second,  the  appropri- 

|  J  *  J 

ate  magnetic  field  Bft..  These  problems  are  both  implicit 

requiring  the  solution  of  a  nonlinear  system  of  equations 

arising  from  the  integration  filters  selected  for  space  and 

» 

time.  These  problems  are  also  both  coupled,  because  BQ« 

«  5  BJ 

depends  on  J  (E_,  8„,  VT  )  and  the  mesh  itself.  As  the 

C  Z  j  Z  P  6 

only  area  of  choice  in  the  processing  sequence,  one  has 

» 

therefore  the  option  of  nesting  the  iterations  for  V  .  and 
»  r  J 

cB9j  in  one  of  two  ways:  velocity  loop  Inside  magnetic 

field  loop,  or  velocity  loop  outside  magnetic  field  loop. 

Since  the  velocity  loop  is  the  most  critical  the  present 

implementation  puts  it  outermost ,  the  latter  option.  This 

has  several  advantages  for  a  subcycled  architecture.  First, 

» 

the  most  recent  solution  of  the  mesh  problem  {  Vp. }  at  a 
previous  call  becomes  a  good  seed  guess  for  the  new  source 
configuration  as  well,  since  it  combines  all  the  information 
about  the  system.  In  contrast,  because  the  mesh  is  very 
sensitive  to  the  local  JxB  stresses,  the  new  velocity  solu¬ 
tion  can  be  far  the  previous  one  in  the  case  of  arbitrary 
_J_,  which  arises  when  the  velocity  loop  is  inside  the  mag- 

C  Z  J  “ 

netlc  field  loop.  Second  the  Jacobian  of  the  mesh  problem 
is  more  nearly  band-diagonal  if  the  magnetic  field  is  assum¬ 
ed  given.  In  contrast,  the  Jacobian  for  the  magnetic  prob¬ 
lem  is  lower-Hessenberg  for  a  given  velocity  field.  Third, 
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the  iteration  of  the  mesh  problem  for  fixed  {  Tt.,  8  .  , 

E  ,}  requires  a  fresh  set  of  densities  n_.,  intermediate 
C  ZJ  ,  C  Ij  » 

velocities  cSj»  thermal  gradients,  ionization  states  cPj » 

and  all  active  transport  coefficients  at  each  evaluation  of 
the  acceleration.  In  contrast,  the  iteration  of  the 
magnetic  problem  requires  fewer  fresh  transport  coeffi¬ 
cients,  constant  densities  and  thermal  gradients,  and  no 
equation  of  state  updates.  Fourth,  the  mesh  problem  has 
more  latitude  in  the  boundary  conditions  -  momentum  flux 
being  legislated  only  if  the  Lagrangian  zone  intercepts  3ome 
fixed  wall  or  the  origin  in  radius.  In  contrast,  the  mag¬ 
netic  field  problem  has  a  firm  boundary  condition  at  the 
origin  which  determines  the  preferential  direction  of  search 
as  the  integral  equation  is  resolved. 

In  view  of  the  considerations  just  set  forth,  the  compu¬ 
tational  advantages  of  putting  the  velocity  loop  outermost 
are  : 


(1)  fewer  total  iterations,  because  the  mesh  problem 
can  be  solved  by  "corrector  only"  methods  and  the  magnetic 
problem  resolves  quickly  for  a  wide  range  of  "seed"  fields 


cJz  j  ’ 


(2)  more  accurate  applicable  matrix-based  methods  for 


V  .  because  the  banded  Jacobians  can  be  factored  more  ac- 
r  j 


cur ately , 

(3)  fewer  total  calculations,  because  the  intermediate 
fields  needed  are  fewer  for  the  magnetic  loop,  and 

(4)  more  efficient  root  searching  because  of  the  pre¬ 
ferred  direction  of  solution  in  the  innermost  loop. 


These  advantages  translate  into  a  fully  self-consistent 
solution  of  the  velocity  and  magnetic  problems  in  an  average 


of  5  outer  and  6  inner  iterations  for  a  total  of  30  evalua- 

*  * 

tions  of  the  derivatives  D  V  ,/Dt  and  J  .,  respectively. 

r  J  C  Z  J 

The  magnitude  of  relative  error  in  the  velocity  corrector 

•  8 

solution  is  usually  bounded  above  by  10  ,  over  all  cells, 

while  that  in  the  magnetic  field  corrector  solution  is 
bounded  above  by  10  nearly  full  machine  precision.  The 

errors  in  the  solution  of  these  nonlinear  difference  equa¬ 
tions  can  be  driven  to  underflow  if  desired,  of  course,  but 
the  experience  to  date  is  that  this  level  of  convergence  is 
adequate  until  the  hydrodynamic  timesteps  can  be  lengthened. 

Once  the  mesh  and  magnetic  field  problems  are  resolved 
to  a  preselected  error  bound,  the  remaining  tasks  are  readi¬ 
ly  performed  in  the  order  dictated  by  the  fluid  equations. 
Working  from  the  highest  velocity  moment  downward,  and  from 
the  interior  to  exterior  in  position  domains,  the  first 
calculation  is  that  for  thermal  sources  and  sinks  (7  •  V  , 

P 
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This  is  followed  by  the  calculation  of 


dissipative  transfers  using  the  flux-limited  forms  of  the 

thermoelectric  and  conduction  processes  recently  developed. 

* .  « . 

From  the  foregoing,  the  material  derivatives  T_  and  0  . 


can  be  formed, 
needed)  of  D/Dt 


o,  - - C*IJ  c~ej 

The  next  calculation  is  an  evaluation  (if 
» 

V  from  the  mesh  configuration  and 


appropriate  sources;  usually  it  is  sufficient  to  accept  the 

» 

last  value  evaluated  in  the  root  search  for  the  V  . 

rJ 

solution.  Finally,  the  calculation  of  D  *E  ./Dt  is  done, 

c  z  j 

using  the  acceleration  field  as  a  source  term  and  coupling 

self-cons  latently  to  the  external  circuit  relation  with  a 

-1  6 

relative  error  of  10  .  The  final  output  is  a  very  accu¬ 

rate  picture  of  these  material  derivatives. 
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THE  SPATIAL  STRUCTURE  OF  CHANGE 


As  an  example  of  the  kind  of  processes  one  must  resolve 
within  the  plasma  load.  Figures  II. 2,  II. 3*  and  II. 4  have 
been  extracted  from  the  early  phases  of  an  implosion  calcu¬ 
lation.  All  three  are  depicted  on  the  same  spatial  scale  a 
2.48  ns  after  the  start  of  test  calculation.  While  the 
initial  configuration  was  a  strictly  isothermal  (25  eV) 

Gaussian  density  and  isoelectric  (E  -3  StV/cm)  profile 

* 

expanding  about  r-0.55  at  a  low  velocity  (  |  Br | <  5.0x10  ), 

the  later  fluid  state  intercepted  these  graphs  depicts  a 

velocity  field  (Br  -  Figure  II. 1)  negative  over  most  of  the 

domains  and  a  strongly  skewed  E  profile;  with  the  outer 

z 

plasma  field  of  44.0  StV/cm  on  order  of  magnitude  than  orig¬ 
inally  given.  Owing  to  the  electromagnetic  stresses,  a  very 
slight  crushing  is  distorting  the  outer  domains  of  the  orig¬ 
inal  Gaussian,  while  the  interior  domains  begin  to  execute  a 
rarefaction  fan  in  ignorance  of  the  piston.  The  exterior 
velocities  have  come  down  considerably  and  the  electron 
temperature  on  the  surface  is  beginning  to  rise. 

2 

From  this  basic  configuration  {n_,  B  .  E_ ,  and  Tt}, 

X  i  Z  6  1 

Figure  II. 2  is  completed  with  the  self  consistent  current 
density  profile  as  resolved  from  the  integral  equation  dis¬ 
cussed  at  length  above.  This  structure  is  essentially  bar¬ 
ren  of  any  thermoe lectr i c  effects  and  shows  the  early  forma¬ 
tion  of  a  peak  in  J  ,  drawn  in  the  domain  where  E  decays, 
due  to  magnetic  insulation  effects.  The  peak  value  is  typi¬ 
cal  of  the  current  density  one  might  expect  in  a  stat ionary 
plasma  for  these  conditions  and  is  contrasted  clearly  with 


2.  T.(r)  is  not  shown  in  Figure  II. 2,  being  essentially 

indistinguishable  from  T@(r). 


that  current  actually  drawn  as  one  rinds  more  magnetic  insu¬ 
lation  in  the  exterior.  Note  also  in  this  context  coinci¬ 
dence  of  the  peak  in  J  with  the  sharp  negative  excursion  in 

Z 

Bp  at  r-0.67  cm.  The  relative  magnitude  of  current  densi¬ 
ties  interior  to  the  main  peak,  and  the  more  gentle  downturn 
in  magnitude  over  the  interval  [0.M9,  0.61]  in  r,  reflects 
the  same  magnetic  insulation  process  obtained  for  a  gentler 

velocity  field  gradient  in  a  smaller  interior  |E  |. 

z 

The  next  graph  [Figure  II. 3]  shows  the  spatial  structure 
of  the  important  source  terms  in  the  electron  temperature 
equation.  The  radiative  loss  profile  Qpad,  bein8  dominated 
by  thin  recombination  continiuum  and  Bremmstrahlung  and 
hence  reflecting  the  density,  is  broad  in  space.  It  also 
proclaims  its  thin  character  by  exhibiting  a  mild  kink  where 
the  density  compression  is  beginning.  In  contrast,  the 
ohmic  heating  profile  Q .  _  is  clearly  peaked  at  the  inner- 

J  •  £ 

most  "insulation  front"  pointed  out  above;  while,  at  the 
exterior,  the  corresponding  insulation  front  peak  has  blend¬ 
ed  with  the  general  upturn  in  ohmic  heating  associated  with 

the  high  E  ,  low  density  region,  the  so-called  corona  plas- 
z 

ma .  At  this  point,  therefore,  the  electromechanical  shock 
has  just  begun  to  emerge  from  the  "vacuum"  and  to  accelerate 
the  load.  The  structure  of  the  shock  heating  profile  Q 

pd  V 

clearly  exhibits  this,  underwriting  as  it  were,  the  wider 
structure  in  the  ohmic  heating  curve.  The  acceleration 
profile  | DV  /Dt |  completes  the  graph  and  shows  a  minimum 
coincident  with  the  slowest  portion  of  the  load.  At  this 
minimum  the  contest  between  Incoming  E  with  its  attendant 

Z 

J  B.  stress)  and  the  pressure  gradient  is  played  out;  while 
z  0 

inside  it,  the  rarefaction  fan  dominates  the  motion. 


f  /  / 


When  all  these  effects  are  Included  with  the  other  dis¬ 
sipative  processes  of  thermal  and  electric  diffusion  the 
plasma  response  contained  is  that  of  Figure  II. 4.  The  re¬ 
maining  material  derivatives  are  shown  and  one  can  see  read¬ 
ily  the  various  regions  of  change  and  the  primary  sources  in 
each.  Outside  the  slowest  velocity  (r>0.6),  the  ion  and 
electron  temperatures  are  being  driven  apart  by  ohmic  heat¬ 
ing,  but  both  are  rising  due  to  compression  and  ohmic  ef¬ 
fects.  In  the  dense  interior,  radiation  and  expansion  are 
cooling  the  fluid,  making  field  diffusion  marginally  faster; 
but  in  the  underdense  interior  a  slight  preponderance  of 
ohmic  heating  arising  from  the  weak  initial  field  begins  to 
warm  the  plasma.  The  material  derivative  of  E_,  on  the 
other  hand,  shows  the  processes  of  field  penetration.  Just 
inside  the  peak  in  the  primary  source  of  incoming  E_  is 

the  term  —  (  B)Bol  which  represents  an  inductive  recoil  — 

T  a 

allowing  more  Ez  to  the  conductor  when  the  flow  field  rises 
to  produce  more  Insulation.  This  recoil  effect  is  superpos¬ 
ed  with  normal  diffusion  and  in  fact  dominates  the  peak  E_ 

11  Z 

shown  at  r-0.65.  On  the  exterior,  where  the  gradient  8 

r  r 

and  the  gross  acceleration  are  somewhat  weaker,  the  diffu¬ 
sive  term  proper  is  competitive  as  is  al30  the  cancellation 

A 

rate  term  E*  fcnE;  but  at  the  shock  itself  the  E  is  moving 

T  Z 

past  the  moving  fluid  predominantly  because  of  the  ac¬ 
celerations  it  gives  the  fluid.  As  the  velocity  field  is 
boosted  inward,  it  would  appear  that  the  E  profile  will 

Z 

simply  follow  this  acceleration  front  along  and  completely 
penetrate  the  load. 

In  the  calculation  of  this  structure  of  change,  it  is 

obvious  that  an  accurate  resolution  of  J_  is  central  to  all 

z 

of  the  important  physics  Just  discussed:  the  insulation 


front,  the  recoil  source  for  E  ,  and  the  regions  of  ohmic 

z 

heating.  It  is  very  important  to  know  the  interplay  of  all 
the  effects  discussed  above  and  to  incorporate  the  non-local 
character  of  the  magnetic  insulation  process  in  the  field 
diffusion  equation. 


CHAPTER  III 


THE  INTEGRATION  OF  THE  FLUID  MESH  TRAJECTORIES 


In  the  original  formulation  of  this  implosion  model, 
correction  of  the  mesh  was  intended  to  occur  outside  the 
subcycle  integration,  allowing  the  explicit  calculation  of 
mesh  positions  in3ide  the  subcycle  derivative  routine  itself 
for  any  input  time  argument.  However,  this  scheme  allowed 
too  much  noise  to  enter  the  velocity  profile  because  the 
feedback  from  thermal  to  flow  reservoirs  tended  to  come  too 
late.  As  a  means  of  controlling  such  noise  one  is  naturally 
led  to  a  predictor /corrector  scheme,  of  necessity  fully 
implicit.  On  the  other  hand,  within  a  3ubcycle  process,  the 
totally  arbitrary  time  coordinate  voids  the  use  of  customary 
formulas  for  the  mesh  predictor /corrector  because  the  time 
levels  are  no  longer  evenly  spaced. 

A.  CONTINUOUSLY  VARIABLE  TIMESTEPS  AND  FILTER  THEORY 


In  order  to  address  the  problem  of  uneven  meshes,  con- 

i- 1  i  * 

sider  a  sequence  of  time  levels  (  t<  t<  t)  together  with 

some  corresponding  velocity  and  acceleration  values  {  1  1u, 

i  *  i-1  i  *  ,  * 

u,  u,  a,  a,  a}.  The  forward  most  time,  t  ,  can  be 

parameterized  by  a  step  size  ratio 


e  =  (  t 


1t)/(it  -  i_1t) 


relative  to  the  step  fixed,  defined  in  the  history  array,  h 


If  one  seeks  an  integrating  filter  of  the  form 

«u  -  aJu+A*-1u  ♦  h(  eB_1  *a  +  BQa>B}_1  a)  ,  [III .  1  D 


10U 


how  can  the  five  coefficients  be  related  to  the  error  char¬ 


acteristics  and  the  timestep  ratio  e?  It  is  useful  to  ad- 
dreis  this  topic  with  a  straightforward  extension  of  the 
transfer  function  as  commonly  defined  in  Fourier  or  z-trans- 
form  theory.  Let  the  accelerations  be  limited  to  harmonic 
components  a- 1 AT  e^ u ( t- * t )  and  seek  a  particular  solution 
u- 1 AQut eJ “ ( t- 1 1 )  ,  the  filter  is  then  charcterized  uniquely 
by  the  ratio 


Aout  .  •  ee,  <•  B0  * 

A .  z  .  - 1 

in  z  -  Aq  -  A i z 


[III  .2] 


where  Z  =  e^ul1. 

For  a  perfect  integrator  [III. 2]  must  equal  1/jw  and  the 
selection  of  coefficients  {A0,  A,,  B- , ,  B0,  B,}  must  be  a 
compromise  process  which  produces  optimal  accuracy  and 
stability  as  a  function  of  the  step  ratio  e.  Let  the  phase 
g=uih  be  an  expansion  parameter  and  require  that  the  perfect 
integrator  be  approximated  through  some  order  in  g .  One 
must  have 


1 


j  ? 


*  E(0[gn]) 


[III  .3] 


where  the  error  term  E  is  bounded  by  gn,  n£3  if  one  wishes  a 

r  i  r  £ 

second  order  or  higher  integrator.  Since  z  -e  the 
sequential  orders  in  g  arise  by  expanding  the  right  hand 
side  of  III. 3  and  equating  real  and  imaginary  parts.  This 
yields  at  each  order  a  new  element  in  the  set  of  constraints 
on  the  coefficelnts  in  the  formula,  viz. 


A1  ♦  e  (1  -  B_1  ) 


BQ  -  B1  -  e2  -  2B_1  +  e ( 1  -  B_, ) 

(2e  ♦  3  e  3  -  e A ^  )  ( 1  *  )  -  (1  -e)A2 

B- 1  "  6e ( c  +  rrn  ♦  ea1  )  * 

These  expressions  already  contain  some  eliminations, 
e.g.,  A 0 ,  but  as  a  set  they  specify  all  of  the  formula  as  a 
function  of  kl.  In  the  case  that  e-1  these  constraints 
reduce  to  a  well-known  class  of  formulas  due  to  Hamming. 
The  new  extension  here  is  to  arrange  the  theory  in  such  a 
way  as  to  produce  uniform  truncation  error  (oCt3])  in  the 
integrating  filter  for  a  continuum  of  timestep  ratios  e. 

As  in  the  fixed-step  case,  the  choice  of  A,  represents  a 
trade-off  between  accuracy,  stability  and  noise  propagation. 
The  introduction  of  the  c*  constraint  could  provide  a  unique 
function  Aj(e)  but,  after  some  investigation  of  the  quartic 
equation  involved,  it  appears  the  real  solutions  with  the 
property  At(0)>0  may  not  exist.  This  property  is  importrant 
because  it  assures  precise  continuity  of  the  velocity  value 
*u  forward  from  iu,  i.e.,  the  filter  collapses  to  an  identi¬ 
ty  transformation  as  e— »0  when  Ax  also  vanishes.  A  second 
point  in  favor  of  avoiding  the  c*  constraint  is  that  such 
high  accuracy  usually  pushes  one  onto  the  stability  bound¬ 
ary.  In  examining  the  present  fomulas  over  the  entire  A,/e 
domain  for  those  cases  with  At-±1  and  e-1,  the  response 
curve  showed  a  small  imaginary  part,  corresponding  to  a 
better  approximation  to  the  perfect  integrator.  These  re¬ 
gions  are  known  stability  boundaries  for  the  fixed  step 
formulas  of  Hamming,  therefore,  such  high  accuracy  can  be 
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traded  here  in  order  to  guarantee  stability.  A  third  point 
in  favor  of  leaving  A,  free  is  the  possible  control  of  noise 
propagation.  A  priori  it  is  of  interest  to  produce  an  inte¬ 
grating  filter  with  a  sharp  or  even  tunable  low  pass  re¬ 
sponse  because  the  local  hydro  cell  sizes  determine  those 
frequencies  too  high  to  be  physical.  The  choice  A^e)  may 
have  some  bearing  on  this  question  as  well.  The  character 
of  the  response  curves  will  be  discussed  at  more  length  in 
the  next  section. 

The  filters  compared  in  later  sections  differ  in  only 
one  respect  -  the  functional  dependence  A  ,(e).  They  are 
arbitrarily  designated  Class  A  and  Class  B,  representing  two 
qualitatively  distinct  types  of  integrating  filter. 

The  Class  A  filter  is  generated  by  the  rule  A  t  -  0.25  e 
and  is  characterized  on  (0<e<1.2)  by  three  features: 

(1)  coefficients  of  largest  magnitude  remaining  posi¬ 
tive  for  all  e, 

(2)  eB_i  uniformly  smaller  than  B0  for  all  e, 

(3)  B  t  >  0  over  most  of  the  e  domain. 

*  i 

This  particular  filter  does  not  mix  a  and  a  with  equal 
weights  until  e>1.2  and  has  ia  predominating  by  a  factor  of 
2  or  more  for  e< 1 . 

The  Class  B  filter  is  generated  by  a  more  complicated 
rule : 

A i  -  [-0.4  +  0.2  f ( e  )  ]  e 

with  f(e)  a  smooth  step  function  [0<e<1.2]  bounded  above  by 
1  and  below  by  0.  This  rule  is  selected  to  track,  roughly, 
a  peak  in  the  value  of  eB_i  at  small  e  until  the  maximum 
blends  into  a  general  rise  for  e  -  0.45.  Thereafter  the 
rule  keeps  A,  well  above  the  stability  boundary  at  -1  for 


e< 1.2.  In  comparison  this  Class  B  filter  is  charcterized  on 
0<e<l .  2  by 

(1)  positive  eB_1  for  all  e, 

(2)  eB,  exceeding  B„  for  larger  e>1.05, 

(3)  B,<0  over  the  domain  studied. 

Hence  the  Class  B  filter  tends  to  use  the  history  as  a 
derivative  estimate  for  the  surge  (B„>0,  B^O),  while  the 
Class  A  filter  tends  to  average  previous  accelerations.  The 
following  sections  are  devoted  to  the  comparison  of  these 
integration  filters  and  to  the  implications  of  these  results 
for  the  hydrodynamic  code. 

B.  RESPONSE  STABILITY  AND  ACCURACY 

If  the  transfer  function  is  decomposed  into  real  and 
imaginary  parts  one  obtains  a  picture  of  the  propagation 
characteristics  of  the  filter  on  the  frequency  domain  defin¬ 
ing  the  time  dependence  of  the  acceleration ,  a(t).  The  more 
useful  object  however  is  the  scaled  difference  of  the  trans¬ 
fer  function  from  the  perfect  integrator  1/Ju»,  i.e.,  H(c)  = 
JCAout/Ain  •  A  comparison  of  Rg  H(c)  and  iffl  H(c)  for  vari¬ 
ous  values  of  e  (between  the  two  classes  of  integrators 
defined  above)  is  a  helpful  summary  of  their  differences. 
The  typical  frequency  response  one  finds  is  illustrated  in 
Figures  111.(1-2)  for  the  Class  A  integrator  and  in  Figures 
III .  (3-1*)  for  the  Class  B  integrator.  The  plots  show  the 
real  part  (R)  imaginary  part  (I)  and  modulus  (M)  of  H(c)  at 
two  values  of  e  in  each  case.  For  e-1.1  the  Class  A  formula 
has  a  slightly  broader  pass  band  in  the  real  part  but  a 
flatter  low  frequency  response.  The  Class  A  formula  shows  a 
larger  magnitude  in  the  imaginary  part  for  all  but  this 
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Figure  III.l 
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Figure  III. 3 


imaginary  part  is  negative  at  high  frequencies.  For  e-0.5 
the  Class  B  formula  retains  its  rather  slowly  decaying  low 
pass  character,  while  the  Class  A  formula  exhibits  a  mild 
rise  in  the  real  part  before  assuming  the  usual  downturn  at 
higher  frequencies.  Both  integration  filters  have  similar 
negative  imaginary  parts.  From  this  comparison  (which  per¬ 
sists  for  all  e)  one  might  expect  the  Class  A  formula  to  be 
slightly  more  susceptible  to  noise  but  more  accurate  and  the 
Class  B  formula  to  be  more  likely  to  propagate  growing  para¬ 
sitic  solutions  due  to  the  more  positive  imaginary  parts. 

A  second  analytic  tool  is  the  relative  stability  criter¬ 
ia  for  each  of  the  integrators.  This  test  is  somewhat  limi¬ 
ted  in  its  scope  because  the  root  conditions  apply  only  to 
even  meshes  and,  while  solutions  for  ei<1  exist,  the  inter¬ 
pretation  is  not  particularly  clear.  The  first  root  condi¬ 
tion  is  independent  of  timestep  and  involves  only  A,  and  A0. 
It  requires  essentially  that  the  denominator  of  the  transfer 
function  possess  no  zeros  for  z>1.  For  the  Class  A  formula 
the  roots  at  e-1  are  z-{1.0,  -0.25);  for  Class  B,  z-{1.0, 
0.2}.  In  both  cases  absolute  stability  is  assured.  The 
second  root  condition  involves  the  phase  also  and  speaks  to 
the  stability  of  propagated  errors,  the  relative  stability 
of  the  scheme.  Since  phase  values  greater  than  0.2  are  not 
likely  to  be  used  on  a  single  step  c*0.2  is  taken  as  the 
phase  in  the  test.  The  corresponding  roots  at  e-1  for  the 
Class  A  formula  are  z-{1.219,  -0.2256};  while  for  class  B, 
z«{1.221,  0.2091}.  Again  in  both  cases,  relative  stability 
seems  assured  with  Class  A  perhaps  a  bit  more  stable  as 
expected . 

It  is  a  well-known  result  that  stable  consistent  methods 
are  convergent  integrators  for  even  timesteps.  Certainly 


one  expects  this  property  to  generalize  smoothly  for  the 
continuously  variable  step  methods  examined  here.  The  best 
kind  of  test  is  therefore  an  integration  of  some  known  and 
(preferably)  relevant  problem.  For  the  integration  of  im¬ 
ploding  trajectories  the  natural  test  bed  is  the  Gaussian 
implosion  a  self-similar  isothermal  (in  space)  fully  time 
dependent  hydrodynamic  flow.  To  summarize  the  fluid  densi¬ 
ty,  temperature,  Lagrangian  position  and  flow  velocity,  and 
acceleration  are  specified  by: 


2  ,  2, 


n  ( r , 

t) 

"0 

e 

s 

\  V  / 

irr 

» 

T  ( t ) 

- 

Vu( 

T  )  , 

r  ( t ) 

m 

r  ( 0 ) 

1/2,  . 

U  (  T  ) 

( f  or 

rg(t)  also), 

V(t) 

m 

V(0) 

u‘1/2 

(x) 

(at 

r  ( t ) ) , 

V(tj 

m 

a  ( 0 ) 

u~3/2 

(t) 

( at 

r  ( t ) )  , 

[III. 4] 


where  the  generating  function 

Y  +  [1  +  (1  +  Y  )t]2 

u(t)  -  — - 7 - r-1 - - -  [III. 5] 

E 

depends  on  the  initial  velocity  scale  and  kinetic  to  thermal 
energy  ratio, 


K  ^  t  /  2 


8  T0/»r^ 


Kj  -  V(0)  ♦  2/rg(0) 
a  (  0  )  -  2  Tn /mr  ( 0  )  . 

0  3 


[III  .6] 
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The  test  consists  in  selecting  K^,  Y£  and,  using  the 
filter  [  1 1 1 . 1  ]  with  V(t)  calculated  analytically  from 
[III. 4-6],  comparing  the  known  velocity  to  the  filter  output 
and  the  known  position  to  the  corresponding  second  intergal 
of  [III.1 ],  viz. 


*r  -  *r  +  eh(  *u  +  *u) /2 


[III. 7] 


This  kind  of  check  tests  only  the  filter,  the  analytical 
accelerations  remove  any  errors  propagated  when  the  filter 
is  combined  with  some  algorithm  for  corrector  convergence 
e.g.,  regula-f alsi .  A  fully  implicit  scheme  for  this  prob¬ 
lem  does  exist  and  will  be  examined  later  if  it  is  of  inter¬ 
est.  As  a  second  exercise  the  acceleration  is  perturbed  by 
a  sinusoidal  ripple,  a«aflu(t)  +  isinwt.  This  allows  a 
check  of  the  frequency  response  in  terms  of  complete 
throughput  as  opposed  to  the  single-step  analysis  above. 

The  timestepping  for  these  checks  is  controlled  by  a 
relative  tolerance  for  single  step  velocity  changes  which 
progressively  slows  the  step  as  the  turning  point  is  reached 
due  to  ever-increasing  accelerations.  For  very  small  toler¬ 
ances  the  step  usually  sinks  to  a  given  floor  sufficient  to 
limit  the  number  of  steps  over  a  full  implosion  to  about 
30,000.  The  relative  error  in  r  and  V  is  reported  every  0.1 
ns  —  y£  and  ^  being  appropriate  to  60  ns  overall  implo¬ 
sion/expansion  times.  The  time  supremum  (in  absolute  value) 
of  r  and  V  errors  at  full  step,  and  1 / 3  ~ »  1 / 2  —  ,  2/3-step 
test  subcycles  is  also  accumulated  on  every  push. 

The  results  are  summarized  in  Tables  III.1  and  III. 2. 
The  Class  A  formula  is  decidedly  superior  on  all  counts  over 
a  wide  domain  in  the  implosion  phase  space  ( Y^ ,  K ^ ) .  The 
Class  B  formula  wins  only  in  the  suppression  of  very  high 


TABLE  III.l 


Gaussian  Implosion  Trajectory  Integration 
at  Single-step  6V/V  =  5.0 ‘10“ 4 

6x  -  xfj^ter  Xanalytical 
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,(3V,S 
sup  (— ) 
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5. 60*10-2 
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9. 90  *10“ 2 
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0.121 

(-0.1,  0.01) 

1.01 *10" 4 

0.304 

1 . 55  *10-4 

0.553 

(-0.1,  0.001) 

7.59*10-4 

0.558 

1.11*10~3 

1.137 

§ 

Relative  velocity  errors  are  dominated  by  the  behavior  near  the  turning 
point,  elsewhere  they  are  comparable  to  ( 6r/r) . 


TABLE  III. 2 


Noisy  Gaussian  Implosion  Trajectory  Integration 
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frequency  noise  above  the  Nyquist  frequency  and  hence  of 
only  minimal  interest.  The  Class  A  formula  exhibits  isolat¬ 
ed  resonances  in  this  domain  while  the  Class  B  first  damps 
them  and  then  admits  a  single  large  one. 


CHAPTER 


SNOWPLOW  FORMATION  - 

At  the  close  of  Chapter  II 
material  derivatives  of  temp 
pointed  out  some  of  the  funda 


density  peak  r  .  ,  while  in  the  last  configuration  it  has 

n  |  p  6  d  K 

moved  well  inside  this  peak.  The  slowest-veloc ity  point 
r6max'  and  the  half “amp  1 i tude  velocity  point  (relative  to 

the  outer  surface  velocity)  r.  _  are  also  seen  to  move  in  a 

similar  manner,  overtaking  the  calculation  zones  as  part  of 
a  general  acceleration  or  "boost”  front  which  is  propagating 
through  the  plasma. 

A  more  detailed  graph  of  the  model  plasma  at  18.8  n3 
(Figure  IV. 2)  shows  the  fully  developed  snowplow  configura¬ 
tion.  The  current  density  now  peaks  more  sharply,  at  a 
deeper  radius  and  the  width  of  the  current  carrying  layer 
has  doubled  over  that  in  the  intial  state.  The  progress  of 
the  boost  front  is  marked  by  the  slowest  velocity  point 
(0.476  cm)  as  compared  to  0.62  cm  at  2.48  ns.  In  step  with 
this  are  both  and  r1/2B  at  ~°*52  cm  (in  from  -0.66  cm 

at  2.48  ns)  and  the  velocity  profile  for  r>0.5245  cm  is 
flattened  in  conjunction  with  a  quasi-equilibrium  configura¬ 
tion  of  the  stresses.  The  ion  and  electron  temperatures 
show  a  peak  at  the  location  of  the  compression  max'..aui  which 
is  in  turn  essentially  coincident  with  the  peak  acceleration 
and  located  just  outside  the  boost  front,  cf.  Figure  IV. 3. 
The  minimum  acceleration  magnitude  lies  inside  the  minimum 
velocity  magnitude  as  in  the  earlier  graphs  and  both  lie 
inside  the  temperature  peak  associated  with  the  shock.  The 
radiative  loss  profile  (cf.  Figure  IV. 3)  now  reflects  this 
in  its  triple  peak:  the  interior  peak  in  T  ,  the  density 

maximum,  and  the  product  of  decaying  density  and  rising 
temperature  in  the  near  corona.  It  is  still  dominated  by 
thin  continuum  and  Br emms trahlung . 

The  field  penetration  process  is  still  dominated  by  the 
recoil  mechanism  -O  B  )Bn  at  the  boost  front  but  is  given 


by  a  mix  of  effects  inside  and  outside  this  point.  The 


configuration  of  J,  E,  and  Bp  has  not  changed  its  character 


since  its  formation  at  early  times,  cf.  Figure  II. 2.  The 
structure  in  the  vicinity  of  the  field  decay  radius  r-,/2E  is 


just  moving  across  the  plasma  load  dragging  the  E_  field 

z 


inward  at  a  rate  fast  compared  to  that  of  pure  diffusion. 
Spatially  the  penetration  can  be  characterized  by  an 
invariant  ordering  of  radii  with  perhaps  time  dependent 
relative  separation  among  them.  Innermost  is  the  minimum  of 


acceleration  magnitude  r 


max  6  ’ 


next  outermost  the  minimum 


flow  speed  r  o  followed  by  the  field  decay  radius  r, 

max  ts  1  /  tb 


and  completed  by  the  onset  of  quasi-equilibrium  flow  r 
The  J  profile  and  the  velocity  half  amplitude  r 
oscillate  around  r 


QSV  * 

1/26  tend  t0 
,  .  so  that  no  one  ordering  of  these 


1  /2E 

radii  is  appropriate. 

The  onset  of  quasi-equilibrium  flow  at  r 


QSV 


can  be  ex¬ 


pected  to  continue  into  the  corona  plasma.  The  linear  equi¬ 


librium  near  r__„  arises  when  V(nT)-oB(E 
Uo  V  z 


♦  Eth  ♦  B  B ) ;  the 


saturated  equilibrium  will  onset  when  V(nT)  -B 


enecs 


obtains  as  Ez  grows  in  the  "vacuum"  region  and  BB  is  held  to 


moderate  levels  by  the  cap  |sj  <  |BqS|  and  the  saturation  of 
J  at  larger  radius.  The  complete  ordering  of  radii, 
apparently  invariant  during  the  early  run-in  or  boost  phases 


13  cnen  lrmax8  < 


<  r_  }  . 

Corona 


The 


max  B  1 /2E  QSV 

calculation  of  the  equilibrium  density  profiles  falling  on 

this  domain  may  prove  a  valuable  modeling  asset  both  in 

simple  applications  and  in  the  smoothing  and  rezoning  of 


implosion  problems.  In  this  connection  however  it  is  worth 
noting  that  the  full  nonlinear  complexity  of  the  thermoelec¬ 
tric  field  must  be  included  as  can  be  seen  in  Figure  IV. 4. 
Any  such  calculation  of  a  quasi-equilibrium  {hj,  B,  J}  given 


typical  {T  ,  6,  E  }  will  find  E  .  a  30-50  percent  effect  as 

6  Z  tn 

one  approaches  the  corona.  Note  also  the  impact  of  the 
thermoelectric  field  on  the  structure  of  E.  The  main  peak 
in  E  is  due  to  recoil  and  diffusion;  the  inward  foot  is 
mainly  due  to  recoil;  but  the  outer  wing  comes  from 
and  the  magnitude  of  the  effect  is  beginning  to  approach  the 
saturated  or  flux-limited  domain  (discussed  in  Section  II. C) 
at  the  outer  radius  of  this  calculation.  One  may  perhaps 
take  some  consolation  in  the  fact  that  the  saturated  domain 
of  Eth  does  not  have  a  great  impact  on  E  but  the  thermo¬ 
electric  effect  is  rather  important  in  these  plasmas. 

The  completion  of  this  rapid  penetration  of  E  into  the 

z 

plasma  load  has  not  yet  been  calculated  but  the  limiting 
factors  will  be  the  disappearance  of  the  sharp  gradient  in 
velocity  and  the  weakening  of  the  accelerations  generally  as 
the  insulation  front  is  allowed  deeper  into  the  load.  If 
what  occurs  here  is  corroborated  by  wider  numerical  exeri- 
ence  then  a  general  picture  of  the  boost  phase  of  a  run-in 
can  be  formulated.  One  would  expect  an  early  field  pene¬ 
tration  determined  by  a  contest  of  rates  and  the  initial 

density  profile  —  E  (with  J  peaked  at  an  insulation 

z  z 

front)  will  penetrate  the  load  until  the  local  velocities  of 
the  interior  are  comparable  to  the  outer  surface.  At  this 
time  the  penetration  will  become  mainly  diffusive  —  at  rate 
characteristic  of  the  radiatively  cooled  interior  plasma  not 
the  hot  corona  however.  When  .the  penetration  settles  down 
into  diffusion,  the  width  of  the  effective  current  carrying 
layer  may  be  quite  large  although  the  stresses  may  not  be 
reduced  proportionately  because  of  the  peaked  nature  of  the 
J  profile  near  the  insulation  front.  The  investigation  of 
this  process  is  therefore  a  high  priority  for  future  work. 


CHAPTER  V 


SUMMARY  OF  PART  II 

In  the  foregoing  discussion,  the  structure  and  perfor¬ 
mance  of  the  electrodif fus  ive  model  of  a  ID  imploding  plasma 
radiator  has  been  detailed.  Some  modifications  to  the  orig¬ 
inal  formulation  have  become  necessary  in  order  to  achieve  a 
self  consistent  picture  of  the  fluid  motion  and  more  modifi¬ 
cations  are  clearly  needed  to  the  integration  package  in 
order  to  obtain  a  useful  performance  with  respect  to  DNA 
objectives . 

On  the  methodological  side,  a  new  integrator  for  the 
fluid  mesh  has  been  developed  and  applied  successfully  to 
known  trajectories.  Certain  discretization  difficulties  in 
the  Braginskil  thermoelectric  field  model  have  been  identi¬ 
fied  in  the  context  of  E-field  diffusion  and  removed  by  the 
development  of  a  simple  flux  limit  on  the  thermal  gradient 
used  as  a  source  of  this  field. 

From  the  physical  standpoint,  the  novel  E-diffusion 
formulation  has  shown  a  sensible  performance  over  short  time 
intervals  and  resolved  some  fundamental  features  of  the  load 
acceleration  processes  in  a  speculative  but  plausible  init¬ 
ial  condition  for  present  DNA  sponsored  experiments.  The 
novel  equation  of  state  involving  generalized  electron  tem¬ 
perature  0  has  been  shown  a  useful  formulation  as  well. 
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